# Automated solver selection for simulation of multiphysics processes in
porous media

Yury Zabegaev
Department of Mathematics
University of Bergen
Bergen, Norway
& Eirik Keilegavlen
Department of Mathematics
University of Bergen
Bergen, Norway
&Einar Iversen
Department of Earth Science
University of Bergen
Bergen, Norway
& Inga Berre
Department of Mathematics
University of Bergen
Bergen, Norway
Corresponding author. Email<EMAIL_ADDRESS>
###### Abstract

Porous media processes involve various physical phenomena such as mechanical
deformation, transport, and fluid flow. Accurate simulations must capture the
strong couplings between these phenomena. Choosing an efficient solver for the
multiphysics problem usually entails the decoupling into subproblems related
to separate physical phenomena. Then, the suitable solvers for each subproblem
and the iteration scheme must be chosen. The wide range of options for the
solver components makes finding the optimum difficult and time-consuming;
moreover, solvers come with numerical parameters that need to be optimized. As
a further complication, the solver performance may depend on the physical
regime of the simulation model, which may vary with time. Switching a solver
with respect to the dominant process can be beneficial, but the threshold of
when to switch solver is unclear and complicated to analyze. We address this
challenge by developing a machine learning framework that automatically
searches for the optimal solver for a given multiphysics simulation setup,
based on statistical data from previously solved problems. For a series of
problems, exemplified by successive time steps in a time-dependent simulation,
the framework updates and improves its decision model online during the
simulation. We describe the solver selection algorithm, present examples of
how the solver selector tunes the solver during the simulation, and show how
it outperforms preselected state-of-the-art solvers for test problem setups.
The examples are based on simulations of poromechanics and simulations of flow
and transport. For the quasi-static linear Biot model, we demonstrate
automated tuning of numerical solver parameters by showing how the L-parameter
of the so-called Fixed-Stress preconditioner can be optimized. Motivated by a
test example where the main heat transfer mechanism changes between convection
and diffusion, we also discuss how the solver selector can dynamically switch
solvers when the dominant physical phenomenon changes with time.

##  1 Introduction

Processes in porous media involve strong multiphysics couplings, such as
between (subsets of) mechanical deformation, fluid flow, and transport. For
such multiphysics simulations, a large part of the computational cost is
associated with solving systems of linear equations. Popular linear solvers
for multiphysics problems apply preconditioners that split the linear system
into subproblems related to separate physical phenomena. This allows for
tailoring of the linear solver, as illustrated below for our model problems of
poromechanics and coupled flow and transport, both of which are highly
relevant for porous media applications.

Solvers for poromechanics have been an active research field for decades [1].
While direct solvers can be applied for small problems, constraints in CPU and
memory necessitate the application of iterative solvers for application-
relevant problem sizes. A common strategy is to apply block preconditioners;
however, a naive block-Jacobi style splitting into subproblems is unstable and
stabilization is needed to ensure efficient preconditioning; see e.g. [2, 3].
This has led to the development of several splitting strategies [4, 5], and
further to the need to tune parameters within the splitting scheme [6, 7].
Finally, the block preconditioner could be accelerated by an outer Krylov
method.

Similar to poromechanics, coupled flow and transport in porous media involve
interaction between processes with different mathematical structures [8]: The
pressure field that governs the flow is parabolic. In most cases, it is almost
elliptical. Transport can involve components and/or energy and involve one or
multiple fluid phases. It is commonly represented by an advection-diffusion
equation, although one of the transport terms may be ignored depending on the
application. In subsurface engineering, transport is usually dominated by
forced advection, and thus the problem is essentially hyperbolic. Accordingly,
the dominant preconditioning strategy, known as Constrained Pressure Residual
(CPR) [9, 10], involves two stages: First, a pressure equation is derived by a
decoupling technique [11, 12] and solved with an approach that handles the
long-range couplings due to the elliptic components of the solution, commonly
by a multigrid method, e.g. [13, 14]. Second, a local correction is applied,
often in the form of an incomplete lower-upper (ILU) factorization. However,
while CPR has been highly successful for transport problems that are
essentially hyperbolic, it is less effective when diffusion plays an important
role, as will be the case for many thermal problems [15]. In this case, an
alternative approach that acknowledges the parabolic nature of the transport
problem has been shown to be more efficient; see [16] and also [17, 18].
Independent of which strategy is applied, it is again common to accelerate the
preconditioner with an outer Krylov scheme.

From our brief discussion of these two model problems, it is clear that the
design of linear solvers for multiphysics problems entails non-trivial
choices: If a block preconditioner is applied, as is usually the case,
decisions are needed for a splitting strategy, subproblem solvers, and an
outer Krylov method. Each of these involves both categorical choices, i.e.,
which class of method to apply, and tuning of parameters within each solver.
Moreover, as the transport problem indicates, the optimal solver strategy may
change dynamically with the solution. For instance, for flow and transport
problems, changes in well configuration can lead to frequent changes in the
flow field, that may impact solver performance. The strategy should also be
adapted to changes in simulation models, for example varying geometry,
boundary conditions etc. As linear solvers often consume a substantial part of
the overall computational time, these choices are of paramount importance.
While the importance of expert knowledge in selecting solvers relative to the
mathematical structure of the problem is clear, users of simulation tools
often lack direct access to such expertise. Consequently, they have to rely on
heuristics for designing linear solvers.

The aim of this work is to develop an automated framework for solver
selection. It seeks to identify optimal solvers, in terms of computational
time, by integrating expert knowledge where available, with an exploration of
solver options. Moreover, the framework should be able to adapt the solver to
changes in the characteristics of the solution. The framework should be
designed to be seamlessly embedded into the simulation process, incurring
minimal costs in computational time and requiring minimal alterations to the
simulation workflow and software.

We emphasize that we consider the problem of solver selection different, and
downstream, from the problem of development of solvers: Throughout this work,
we assume that for any linear system considered, there is at least one
candidate linear solver strategy, and our task is to make optimal choices
within this strategy. As new solvers are introduced our framework, this opens
up the possibility of computational savings, but it also broadens the range of
options that must be considered in the selection process. Furthermore, the
theoretical study of solvers will often lead to better understanding of
available algorithms, such as permissible ranges for solver parameter values,
that can be used to narrow the choices that must be explored by our framework.

The solver selection problem is a specific application of a more general
algorithm selection problem [19, 20, 21], which also includes a broad range of
other applications such as machine-learning algorithm parameter tuning,
compiler optimization, Boolean satisfiability problems, and search algorithms,
see e.g. [20, 21] for more information. The general task is to find an
appropriate algorithm for a given problem, based on characteristics of the
problem and past performance data for the candidate algorithms. The algorithm
selection methods have found various applications related to numerical
simulations, including, but not limited to, those discussed below.

An established approach to solver selection is formulated as a classification
problem, see, e.g., [22, 23, 24, 25, 26]. Each problem in a training data set
is solved exhaustively by all available solvers, and a machine learning
algorithm learns to map a problem to its optimal solver. This mapping is then
exploited for problems not observed in the training data set. The
classification approach requires that the performance data is available for
all candidate solvers, however, this is not a feasible for problems with
frequently changing dynamics, as illustrated by subsurface transport with
changing well configurations. Incomplete performance data can be dealt with by
viewing the selection process as optimization under uncertainty. A natural
approach is to apply Bayesian optimization, wherein a probabilistic model is
constructed for the target function and used to guide the solver selection.
Works in this direction include [27, 28], which all apply a probabilistic
representation based on Gaussian Processes (GP) [29], or generalizations
thereof.

Enhancements in solver selection can be attained by conceptualizing a
simulation not as a singular process to be optimized, but as a series of
successive problems. Each time step, or each iterative linear solution within
a time step or within a nonlinear simulation, represents a distinct solver
selection problem. This entails a much deeper embedding of the solver
selection process into the numerical simulation, which shortens the feedback
loop and provides frequent information on solver performance. The iterative
nature of time stepping is exploited in [30] to select iterative solvers,
preconditioners, and tune solver parameters. Additionally, the framework can
exploit performance data in databases from previous simulations. Selecting the
linear solver for individual time steps also opens up the possibility of
tuning the solver to changing characteristics of the solution, as was done,
e.g., in [31, 32]. Common to these works is the manual definition of a
problem-specific heuristic for when to switch solvers. Finally, we note that
these characteristics can be learned [33].

Our prime interest herein is to minimize the computational cost of solving
linear system, hence our focus on performance data. Other techniques can also
be applied to obtain information on the system, and thus speed up simulations.
This includes the use of a posteriori error estimates within inexact Newton
schemes, as was done in e.g. [34, 35]. We also note that there has recently
been considerable interest in applying machine learning to tune parameters in
multigrid methods, see e.g., [36, 37, 38].

In this work, we build a solver selection framework which does not rely on a
prior database of solver performance and also adapt to changes in the
solution. To that end, our framework relies on two key ingredients: First, we
consider the selection of a solver for every time step as a separate
optimization problem and update the selection model with newly obtained
performance data. This facilitates automatic experimentation with solver
configurations with a short feedback loop and thus mitigates the lack of
comprehensive solver performance data prior to simulation. Second, we build a
relation between the performance of solver options and the characteristics of
the solution and use this to steer the adaptation between time steps. To the
best of our knowledge, this is the first framework that combines these two
features for the purpose of solver selection. It should be noted that, while
we limit solver selection to individual time-steps in the current work, the
framework is also applicable to solver selection for iterative solutions of
non-linear problems that are not time-dependent.

The rest of the manuscript is structured as follows: Section 2 introduces the
terminology needed to formally define the solver selection problem, while our
methodology is presented in Section 3. In Section 4, we present the governing
equations for the two model problems of poroelasticity and non-isothermal
coupled flow and transport, together with popular linear solver approaches for
these problems. The numerical examples presented in Section 5 show that our
framework is indeed capable of selecting an optimal solver and to adapt the
choice to changing solution characteristics. We also investigate different
approaches to exploring uncertainty caused by scarce data, and show how the
data from one simulation setup can be utilized for another simulation setup to
speed up the selection process for it. Finally, our concluding remarks are
given in Section 6.

##  2 Problem statement

In this section, we introduce the solver selection problem which consists of
three components: The space of available solver configurations, a
characterization of the static and dynamic features of the simulation, and the
reward function.

###  2.1 The solver configuration space

The solver selection framework seeks the optimal solver configuration in a
user-provided solver configuration space which describes all candidate solver
configurations. The structure of this space can be non-trivial, as can be
illustrated by an example for a linear single-physics problem, with a limited
set of solver options: The linear system is solved by either a direct method
or an iterative Krylov solver, taken to be either the Generalized Minimal
Residual (GMRES) or the Conjugate Gradient (CG) method. If a Krylov solver is
applied, it is preconditioned by an Algebraic Multigrid (AMG) method or the
incomplete LU (ILU) factorization. Picking from these options leads to
additional choices to be made, such as the restart parameter in GMRES, the
drop tolerance in ILU and the type of cycle (V or W) in AMG, and so on. The
solver configuration space for this example is visualized in Figure 1(a).

From this example, it is clear that deciding for a specific configuration
involves two principal types of decisions: categorical and numerical.
Categorical decisions entail choosing specific solver components, e.g., the
Krylov scheme and the type of preconditioner, while numerical decisions
involve the setting of real or integer values, e.g., the drop tolerance in
ILU. An important distinction between categorical and numerical decisions is
that the latter has a natural ordering of its values which can be exploited by
a solver selection algorithm. Moreover, there is a finite number of
combinations of categorical choices, while the numerical decisions that
involve real numbers can take infinitely many values. As can be seen in Figure
1(a), the choices have a structure reminiscent of a decision tree, in that a
specific decision may trigger or exclude a later decision: The choice of
preconditioner is only relevant if a Krylov solver was chosen, and the ILU
drop tolerance is not needed in AMG.

We refer by the term solver configuration to the complete algorithm capable of
solving a linear system. The set of all possible combinations of categorical
choices is denoted $\mathcal{A}_{\text{cat}}$, and, for any
$a_{\text{cat}}\in\mathcal{A}_{\text{cat}}$, we let
$\mathcal{A}_{\text{num}}(a_{\text{cat}})$ denote the corresponding space of
numerical choices. The space of solver configurations, $\mathcal{A}$, can be
represented as

$\mathcal{A}=\bigcup_{a_{\text{cat}}\in\mathcal{A}_{\text{cat}}}\mathcal{A}_{\text{num}}(a_{\text{cat}}).$

In our example there are in total seven solver configuration groups, see
Figure 1(b), but only some of the groups contain numerical parameters.

(a)

(b)

Figure 1: (a) Example solver space. The hatched blocks represent categorical
decisions. (b) The seven solver configurations.

###  2.2 The context space

The ideal solver configuration will depend on the properties of the problem at
hand, such as the geometry of the model, the boundary and initial conditions,
the source terms, the model parameters, etc. Knowledge of such factors can be
used to guide solver selection, as has also been shown previously [33]. In
fact, this can be viewed as a mechanism for a user to provide the selection
process with expert insight into the characteristics of a given simulation
model. To that end, a characterization which includes both static and dynamic
features of the problem and solution, can be collected into an array of
numerical values that we will refer to as the simulation context.

We define the simulation context as a space that includes all the possible
values that the simulation context can take,
$\mathcal{C}\subset\mathbb{R}^{m}$, where $m$ is the dimension of the context.
We will provide the content for the specific problems in the corresponding
experiments in Section 5.

###  2.3 Optimization problem

The goal of the solver selection process is to find the optimal solver
configuration $a^{\star}\in\mathcal{A}$ which minimizes the time spent on
solving the linear system and penalizes failed attempts. The optimization
problem is expressed in terms of the reward function
$R:\mathcal{A}\times\mathcal{C}\rightarrow\mathbb{R}$. Given the solver
configuration $a\in\mathcal{A}$ and the context of the simulated model
$c\in\mathcal{C}$, $R$ evaluates the solution efficiency:

$R(a;c):=\begin{cases}-\log T_{\text{sol}},&\text{if successfully solved},\\\
P,&\text{otherwise}\end{cases}$ (1)

where $T_{\text{sol}}$ is the time spent on solving the current problem. The
penalty for a failed attempt, $P$, is set to twice the largest solver runtime
for a given problem, that is:

$P:=-\log\left(2\cdot\max{}T_{\text{sol}}\right)$

The value of the reward function for a given solver configuration $a$ and
context $c$ is not deterministic; therefore, $R$ should be considered a
stochastic quantity. There are in general two sources to the lack of
determinism: First, the context space may not fully describe the features of
the simulation which should inform the solver selection; indeed, if the
context space is void, no such information is available. Thus, for a specific
configuration, $R$ may take multiple values for the same context. Although
this can be alleviated by employing a rich description of the context, such an
approach risks putting a heavy burden on the solver-selection method. Second,
the measurement of the runtime is associated with stochastic measurement noise
caused by other programs executed on the same computer at the time of
measurement. This can in principle be corrected by carefully monitoring the
processing running on the computer, but adjusting for such factors is
technically challenging. A different workaround, to construct the reward
function based on iteration count instead of runtime, is not practical since
candidate solver configurations will have different runtime cost for their
iterations. We chose instead to treat $R$ as a stochastic function and assume
it is described by an unknown stationary distribution. We define the optimal
solver configuration $a^{\star}$ by incorporating the stochastic nature of the
reward function:

$a^{\star}=\operatorname*{\text{argmax}}_{a\in\mathcal{A}}\mathbb{E}\left[R(a;c)\right]$
(2)

where $\mathbb{E}\left[\cdot\right]$ denotes the expectation. That is, we
define the optimal solver as the one which produces the highest reward on
average. In the following, we will denote by sampling the reward function the
process of applying a solver configuration to a problem and evaluating the
corresponding reward.

##  3 Solver selection framework

Since the distribution of the objective function is unknown, we do not have a
functional expression for its expectation. Thus, the formula Eq. (2) defines
the optimal solver, but it does not provide a recipe for how to find it. In
practice, we can only approximate the expectation of $R$ by applying black-box
optimization methods and sampling $R$ with various solver configurations $a$
and contexts $c$. As our interest is problems where each sampling is
considered a computationally costly operation, it is critical that the
optimization method efficiently deals with data scarcity.

Our optimization approach consists of two main ingredients: First, we apply
online learning by assigning, in principle, a new solver for each time step in
a simulation. This provides a short feedback loop for the optimization tool
and allows for rapid adjustment based on newly generated and problem-specific
data from $R$. For simplicity, the solver is fixed for all Newton iterations
within one simulation time step. Second, we use machine learning regression
methods to build a model for predicting $\mathbb{E}\left[R(a;c)\right]$ and
selecting solvers. The solver selection will need to balance exploitation,
where the best solver is selected based on currently available data, with
exploration, which seeks to address lack of data in regions of the solver or
context space.

###  3.1 Online learning

We consider each successive time-step as an independent problem to which we
apply the solver selection algorithm. All linear solutions within a time step
generate performance data that is used to update the solver selection model
during simulation, in what we will refer to as online learning. This allows us
to start a simulation without any performance data, but data from previous
simulations can be included. As the solver selection aims to minimize the time
to solve each linear system, the process will also further the goal of
minimizing the total execution time of the whole simulation.

Compared to a hypothetical setting where the ideal solver is known in advance,
there are two sources of additional computational cost under online learning:
First, with limited performance data available, exploration of the solver
configuration space will necessarily lead to suboptimal solvers being chosen.
This cost can to some degree be alleviated by the selection algorithm, but it
is also an unavoidable consequence of lack of information. Second, time spent
on updating the solver selection model with new performance data and selecting
a solver should be considered part of the overall simulation time. It is
important to note that, while linear system solution time scales with the
system size, the cost of the optimization algorithm will typically scale with
the number of data points available. Thus the cost of optimization will matter
less for large problems where identifying appropriate solvers can be expected
to bring the largest absolute decreases in simulation time.

###  3.2 Machine learning methods

To approximate $\mathbb{E}\left[R(a;c)\right]$ we employ machine learning
regression methods. As discussed in Section 2.1, a solver configuration $a$
consists of categorical decisions of subalgorithms and their corresponding
parameters. Accordingly, we split the solver configuration into
$a:=\left[a_{\text{cat}},a_{\text{num}}\right]$. For brevity, we denote by $x$
a pair of the solver configuration and the simulation context:
$x:=\left[a,c\right]\in\mathcal{A}\times\mathcal{C}$; the reward function Eq.
(1) can take this argument: $R(x):=R(a;c)$. The machine learning algorithms
will rely on performance data, with the set of data points $D$ defined as:

$D=\left\\{x_{i},r_{i}\right\\}_{i=1..n},$

where $r_{i}$ is the value of the reward function sampled for $x_{i}$, and $n$
is the number of data points in $D$. No distinction is made between data from
previous simulations and the current simulation.

For the purpose of improving the accuracy of a machine learning algorithm, the
data set is transformed to have zero mean and unit variance. For numerical
parameters that take real numbers, we split the admissible interval into a
finite set of evenly distributed values and let the optimization choose
between them. In practice, we use $20$ points for each numerical parameter.

Different elements in $\mathcal{A}_{\text{cat}}$ can share components, for
instance both GMRES and CG can use the ILU preconditioner, referring to Figure
1(b). Apparently, this would allow for sharing performance data between solver
configurations and hence a faster reduction in the uncertainty of $R$.
However, while it is easy to measure the time spent in the individual
components of a solver, this information is highly dependent on the
interaction with other solver components and will have limited to no value if
transferred to a different set of categorical choices. Therefore, we disregard
the possibility of such transfer, and instead assign different machine
learning instances to each $a_{\text{cat}}\in\mathcal{A}_{\text{cat}}$,
motivated by [39].

If no data is available for a solver configuration at the start of the
simulation, the first $n_{init}$ data points for each $a_{cat}$ are always
collected via random exploration to provide the machine learning model with
sufficient data to make meaningful predictions. In practice, we have found
$n_{init}=3$ to be a reasonable parameter value.

For the machine learning method trained on data $D$, we denote by
$\mathbb{E}_{D}[R(x)]$ the approximation of the expectation $\mathbb{E}[R(x)]$
defined in Section 2. In the following, we will outline two specific nonlinear
regression methods which differ in the way they treat uncertainty, as the
first incorporates Bayesian statistics, while the second applies heuristics.

###  3.3 Bayesian optimization approach

The classical technique for expensive black-box optimization under data
scarcity is Bayesian optimization [40]. We give a brief overview of it here;
for more information, see [29, 40]. The core idea is to build a probabilistic
approximation of the expected value and variance of the target function, based
on prior data. In practice, the Gaussian Process (GP) statistical model is
applied. The GP algorithm considers a random variable $R(x)$ for each point
$x\in\mathcal{A}\times\mathcal{C}$ and assumes that for a finite number of
points $x$, the values $R(x)$ follow a multivariate normal distribution. The
variance at a point $x$ given data $D$ is denoted $\sigma^{2}_{D}[R(x)]$ and
will increase as we move away from available data, thus the GP representation
allows us to explicitly estimate our confidence in the current model.

The decision-making process of the GP algorithm is guided by an acquisition
function. We utilize the Upper Confidence Bound (UCB) acquisition function:

$\text{UCB}(x)=\mathbb{E}_{D}[R(x)]+\alpha\cdot\sigma^{2}_{D}[R(x)],$ (3)

where $\alpha$ is a parameter; we note that other acquisition functions are
also used, for example, in [28]. The next solver to be evaluated
$\tilde{a}^{\star}$ is selected by maximizing the acquisition function:

$\tilde{a}^{\star}=\operatorname*{\text{argmax}}_{a\in\mathcal{A}}\text{UCB}(x).$

The chosen configuration $\tilde{a}^{\star}$ attempts to approach the true
optimum $a^{\star}$ of Eq. (2) by balancing exploration and exploitation: For
low uncertainty, the first term will dominate the acquisition function and the
solver with the highest expected reward is chosen, thus exploiting existing
data. If the second term in Eq. (3) dominates, the selection will explore
parts of the solver configuration space with high uncertainty. The parameter
$\alpha$ determines how much preference to give to exploration. The space of
functions predicted by the GP algorithm can be controlled by the GP kernel,
which in this work is taken as the radial basis function kernel:

$k(x_{i},x_{j})=\exp{\left(-\dfrac{d(x_{i},x_{j})^{2}}{2l^{2}}\right)},$

where $d(\cdot,\cdot)$ is Euclidean distance and $l>0$ is the algorithm
parameter that corresponds to the length scale of the modeled function. This
parameter is learned during fitting the GP algorithm with the data.

As outlined in Section 3.2, we use separate instances of the GP method for
different combinations of categorical decisions. Each instance is trained with
the corresponding subset of data from $D$. Selection of a new solver
configuration entails evaluating all the GPs and picking the configuration
that corresponds to the best prediction of all the GP instances.

GP is a powerful and established algorithm, which is commonly used in the
solver selection domain [27, 28]. Its computational cost can be significant,
though: The cost of updating the GP with new data is cubic with respect to the
number of data points in $D$, and the cost of selecting a new solver is cubic
with respect to the number of solver configurations to be selected from. Hence
using GP for solver selection may introduce a significant overhead, in
particular for simulations that produce many data points or involve relatively
small linear systems. The performance of GP and the associated overhead will
be tested in Section 5.

###  3.4 Heuristic approach

As a computationally cheap alternative to GP, we here consider a heuristic-
based approach to the black-box optimization. This approach does not apply a
representation of variance, thus we cannot explicitly estimate the uncertainty
and utilize the UCB acquisition function. Nevertheless, we can still
approximate the expectation of the reward function by applying any machine
learning regression model.

To balance between exploration and exploitation we apply an epsilon-greedy
heuristic. With probability $\varepsilon\in[0,1]$, we select a random solver
configuration from $\mathcal{A}$. Otherwise, we make a greedy choice by
selecting the configuration $\tilde{a}^{\star}$ that maximizes the current
approximation of the expectation:

$\tilde{a}^{\star}=\operatorname*{\text{argmax}}_{a\in\mathcal{A}}\mathbb{E}_{D}\left[R(x)\right]$

Note that $\tilde{a}^{\star}\neq a^{\star}$ from Eq. (2) because the former
represents our current belief based on the data $D$, and the latter represents
the true optimum, which is unknown. After each exploration, the probability of
exploring next time decreases: $\varepsilon\leftarrow\varepsilon\cdot\gamma$
where $\gamma$ is a parameter.

The expectation of the reward function is approximated by the gradient
boosting regression algorithm [41]. This machine learning model approximates
data with piecewise constant functions and supports non-linear and
discontinuous target functions. Once fitted, the gradient boosting algorithm
cannot be updated with newly generated data, so the model must be refitted
from scratch after each new time step. The cost of this operation is however
insignificant, as will be shown in Section 5.

As discussed in Section 3.2, separate instances of the gradient boosting
algorithm are used to model the reward for different categorical decisions
$a_{\text{cat}}$. The parameters $\varepsilon$, $\gamma$ and $n_{init}$ are
separate for each instance, and if one algorithm instance decides to explore,
it is selected regardless of how high the predicted rewards are from the other
algorithm instances. Thus, the overall probability to explore is
$1-(1-\varepsilon_{1})\cdot(1-\varepsilon_{2})\cdot\ldots\cdot(1-\varepsilon_{k})=1-\Pi_{i=1}^{k}(1-\varepsilon_{i})$
for $k$ instances of the algorithm.

##  4 Model problems

In this Section, we define the two model problems that will be used to test
the solver selection framework: Deformation of a poroelastic medium and non-
isothermal porous media flow of a single-phase fluid. For each problem, we
present the governing equations, discretization, and linear solver algorithms.
The following information is common for both models: The porous medium is
defined in the domain $\Omega\in\mathbb{R}^{d}$, where $d$ is the spatial
dimension. The domain boundary consists of two disjoint subsets:
$\partial\Omega=\Gamma_{N}\cup\Gamma_{D}$ that correspond to boundary
conditions Neumann ($\Gamma_{N}$) and Dirichlet ($\Gamma_{D}$) boundary
conditions. The problem is time-dependent with $t\in[0,t_{\text{end}}]$.

###  4.1 Poromechanical deformation

The widely used mathematical model for poroelasticity is the quasi-static
linear Biot model. The problem is based on mass conservation and force balance
equations (see, e.g., [42]) and can be stated: Find the displacement
$\mathbf{u}$ and the pressure $p$ such that

$\begin{cases}\nabla\cdot\sigma=\mathbf{0}\\\
\rho\dfrac{\partial\phi}{\partial
t}+\rho\nabla\cdot\mathbf{q}=f\end{cases}\quad\text{in
}\Omega\times[0,t_{\text{end}}]$ (4)

Here $\rho$ is the fluid density, $\phi$ is the rock porosity, $\mathbf{q}$ is
the fluid velocity, $f$ is the mass source term and $\sigma$ is the total
stress. We apply Dirichlet and Neumann boundary conditions on the
corresponding parts of the boundary $\Gamma_{N}$ and $\Gamma_{D}$.

#### Constitutive laws

The fluid velocity is governed by Darcy’s law,

$\mathbf{q}=-\dfrac{\mathbf{K}}{\mu}\nabla p,$ (5)

where $\mathbf{K}$ is the permeability tensor, $\mu$ is the fluid dynamic
viscosity and gravitational effects are neglected. The total stress is defined
as

$\mathbf{\sigma}=\mathbf{\sigma^{\prime}}-bp\mathbf{I},$ (6)

where $\sigma^{\prime}$ is the effective stress, $b$ is Biot’s coefficient and
$\mathbf{I}$ is the second-order unit tensor. The effective stress is
connected to the displacement via Hooke’s law (assuming zero reference stress
and displacement) by

$\mathbf{\sigma^{\prime}}=\mathbb{C}:\mathbf{\epsilon},$ (7)

where $\mathbb{C}$ is a fourth-order tensor stiffness and $\epsilon$ is
defined as follows:

$\mathbf{\epsilon}=\nabla^{s}\mathbf{u}:=\dfrac{1}{2}\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right).$
(8)

We consider the isotropic homogeneous rock, so Hooke’s law can be simplified
to

$\sigma^{\prime}=2\mu_{r}\epsilon(\mathbf{u})+\lambda\nabla\cdot\mathbf{u}\mathbf{I},$
(9)

where $\mu_{r}$ and $\lambda$ are the Lame parameters. In our poromechanics
model, the fluid density and viscosity are constant, while porosity is a
function of pressure and displacement,

$\phi(p,\mathbf{u})=\dfrac{1}{m}p+b\nabla\cdot\mathbf{u},$ (10)

where $\dfrac{1}{m}$ is the specific storage.

#### Discretization

We use fully implicit Euler time discretization. The spatial discretization is
based on cell-centered finite volume schemes. We apply the two-point flux
approximation [43] defined for diffusive scalar problems and the multi-point
stress approximation (MPSA) for vector problems [44]. The spatial
discretization is locally conservative for mass and momentum.

#### Solver configurations

For the given model problem, each time step of the simulation requires solving
a linear system of equations. We solve it with GMRES preconditioned by the
Fixed-Stress approximation of the Schur complement (see [2, 6]). The Fixed-
Stress scheme is briefly reviewed below.

We can rewrite the system (4) in the residual form by moving the source term
to the left-hand side and denoting $r_{\text{mass}}$ and $r_{\text{force}}$ as
the residuals of the mass and force balance equations, respectively. The
corresponding matrix has the following block structure:

$J=\begin{bmatrix}A&B_{1}\\\ B_{2}&C\\\ \end{bmatrix},$ (11)

where $A={\partial r_{\text{force}}}/{\partial\mathbf{u}}$, $B_{1}={\partial
r_{\text{force}}}/{\partial p}$ , $B_{2}={\partial
r_{\text{mass}}}/{\partial\mathbf{u}}$ and $C={\partial
r_{\text{mass}}}/{\partial p}$. The matrix can be factorized by a block $LDU$
decomposition:

$J=\begin{bmatrix}I&\\\ B_{2}A^{-1}&I\\\ \end{bmatrix}\begin{bmatrix}A&\\\
&S_{A}\end{bmatrix}\begin{bmatrix}I&A^{-1}B_{1}\\\ &I\end{bmatrix},$ (12)

where the term $S_{A}$ is the Schur complement with respect to $A$:

$S_{A}=C-B_{2}A^{-1}B_{1}.$ (13)

We approximate $J$ by the block upper-triangular matrix:

$J^{-1}\approx U^{-1}D^{-1}=\begin{bmatrix}A^{-1}&-A^{-1}B_{1}S_{A}^{-1}\\\
&S_{A}^{-1}\end{bmatrix}.$ (14)

The Schur complement should not be computed explicitly as it requires
inverting the matrix $A$, thus an approximation must be sought. The most
common approximation for the poromechanics system (4) is to the so-called
Fixed-Stress preconditioner, which for finite volume discretizations is
defined as follows:

$-B_{2}A^{-1}B_{1}\approx L\dfrac{\rho V}{\Delta t}.$ (15)

Here $V$ is the cell volume and $\Delta t$ is the simulation time step. The
stabilization parameter $L\in\mathbb{R}$ has a significant impact on the
convergence rate of the iterations. In a recent study [7], the analytical
limits for the optimal value of $L$ were derived:
$L_{\text{opt}}\in\left[L_{\text{min}},L_{\text{phys}}\right]$, where

$\displaystyle L_{\text{min}}=\dfrac{b^{2}}{4\mu_{r}+2\lambda}\text{,}\quad
L_{\text{phys}}=\dfrac{b^{2}}{\dfrac{2\mu_{r}}{d}+\lambda}.$ (16)

However, it was shown in [45], that for a given problem, the value of
$L_{\text{opt}}$ depends on the boundary conditions and flow parameters.

###  4.2 Non-isothermal flow in porous media

Assuming single-phase flow in porous media, a model that combines flow and
thermal transport is based on mass and energy conservation equations, which
read, respectively:

$\begin{cases}\phi\dfrac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathbf{q})=f\\\\[10.0pt] \phi\dfrac{\partial h}{\partial
t}+(1-\phi)\dfrac{\partial h_{r}}{\partial
t}+\nabla\cdot(h\mathbf{q})+\nabla\cdot\mathbf{q_{T}}=f_{T}\end{cases}\quad\text{in
}\Omega\times[0,t_{\text{end}}]$ (17)

where $h$ is the specific enthalpy of the fluid, $h_{r}$ is the specific
enthalpy of the rock, $\mathbf{q_{T}}$ is the diffusive heat flux and $f_{T}$
is the energy source. The primary variables are pressure $p$ and temperature
$T$. We apply Dirichlet and Neumann boundary conditions on the corresponding
parts of the boundary $\Gamma_{N}$ and $\Gamma_{D}$.

We use the Peaceman well model [46, 47] to simulate wells; details can be
found in B.

#### Constitutive laws

In the non-isothermal flow model in porous media, porosity is constant and
gravity is neglected. The density is a function of pressure and temperature
while the viscosity is a function of temperature; the details are given in B.
The specific enthalpy of fluid and rock are defined by the following
simplified expressions: $h=\rho c_{v}T$, $h_{r}=\rho_{r}c_{r}T$, where
$\rho_{r}$ is the density of rock and $c_{v}$ and $c_{r}$ are the specific
heat capacities of fluid and rock, respectively. The diffusive heat flux is
defined by Fourier’s law: $\mathbf{q_{T}}=-k_{T}\cdot\nabla T$, where
$k_{T}=\phi k_{Tf}+(1-\phi)k_{Tr}$; $k_{Tf}$ is a thermal conductivity of
fluid, and $k_{Tr}$ is a thermal conductivity of rock.

#### Discretization

We use fully implicit Euler time discretization and finite volume methods for
spatial derivative: Diffusive terms are discretized by two-point flux
approximation, while advective fluxes are discretized by a first-order
upstream scheme.

#### Solver configurations

The non-isothermal flow model is nonlinear and is solved by Newton’s method.
We adjust the time step dynamically so that Newton’s method converges in
around 4 iterations. At each Newton iteration, we solve the full Jacobian
linear system with GMRES. Restricted to each time step, we use as
preconditioner either the Constrained Pressure Residual (CPR) [9, 10] or the
or the Schur complement approximation [16], both of which are briefly outlined
below.

Similar to the poromechanics system from Section 4.1, we rewrite Eq. (17) in
residual form, and define $r_{\text{mass}}$ and $r_{\text{energy}}$ as the
residuals for the mass and energy equations, respectively. Borrowing notation
from [16], the block structure of the Jacobian matrix can be written

$J=\begin{bmatrix}A_{pp}&&A_{pT}\\\ A_{Tp}&&A_{TT}\\\ \end{bmatrix},$ (18)

where $A_{pp}={\partial r_{\text{mass}}}/{\partial p}$, $A_{pT}={\partial
r_{\text{mass}}}/{\partial T}$, $A_{Tp}={\partial r_{\text{energy}}}/{\partial
p}$ and $A_{TT}={\partial r_{\text{energy}}}/{\partial T}$.

The CPR preconditioner is a multiplicative two-step preconditioner. The
algorithm consists of the application of two operators: $M_{1}^{-1}$ and
$M_{2}^{-1}$. For a right-hand side vector $b$, it is defined as follows:

  1. 1.

Precondition using $M_{1}:x_{1}=M_{1}^{-1}\cdot b$

  2. 2.

Compute the new residual: $b_{1}=b-J\cdot x_{1}$

  3. 3.

Precondition using $M_{2}$ and correct: $x=M_{2}^{-1}\cdot b_{1}+x_{1}$

The operator $M_{1}$ approximately inverts the pressure block:

$M_{1}^{-1}\approx\begin{bmatrix}A_{pp}^{-1}&0\\\ 0&0\end{bmatrix},$ (19)

where $A_{pp}^{-1}$ is approximated using one AMG V-cycle. The operator
$M_{2}$ approximately inverts the full system $J$ using an incomplete LU
factorization (ILU). While the CPR preconditioner applied to multi-phase flow
usually includes a decoupling operator, e.g. [11], this was shown not to lead
to additional improvement for the non-isothermal single-phase flow model [16],
and we therefore do not apply any decoupling operator.

The second preconditioner is based on a Schur complement approximation. Again,
we apply the block $LDU$ decomposition to the matrix (18):

$J=\begin{bmatrix}I&0\\\ A_{Tp}A_{pp}^{-1}&I\\\
\end{bmatrix}\begin{bmatrix}A_{pp}&0\\\ 0&S_{T}\\\
\end{bmatrix}\begin{bmatrix}I&A_{pp}^{-1}A_{pT}\\\ 0&I\end{bmatrix},$ (20)

where $S_{T}=A_{TT}-A_{Tp}A_{pp}^{-1}A_{pT}$ is the Schur complement. The
inverse of the block decomposition is given by:

$J^{-1}=\begin{bmatrix}I&-A_{pp}^{-1}A_{pT}\\\
0&I\end{bmatrix}\begin{bmatrix}A_{pp}^{-1}&0\\\ 0&S_{T}^{-1}\\\
\end{bmatrix}\begin{bmatrix}I&0\\\ -A_{Tp}A_{pp}^{-1}&I\end{bmatrix}.$ (21)

Three different preconditioners can be derived by different approximations to
this block inverse: What we will term the full preconditioners include all
three factors, while upper and lower block triangular versions are found by
including, respectively, the two first and the two last matrices. We
approximate the Schur complement as,

$\widetilde{S}_{T}=\phi\dfrac{c_{v}\rho}{\Delta
t}+(1-\phi)\dfrac{\rho_{r}c_{r}}{\Delta
t}+\nabla\cdot(c_{v}\rho\mathbf{q})-\nabla\cdot(k_{T}\nabla)-c_{v}f_{\text{prod}},$
(22)

where $\Delta t$ is the simulation time step and $f_{\text{prod}}$ is the
source term corresponding only to the production wells, see [16] for a
derivation. The preconditioner of Eq. (21) approximates each of the inverses
of the matrices $A_{pp}$ and $\widetilde{S}_{T}$ by one AMG V-cycle iteration.

The quality of the approximation of $\widetilde{S}_{T}^{-1}$ depends on the
ratio between the diffusive and advective terms in Eq. (22), that is, of the
Peclet number: The AMG method can be expected to perform well for low Peclet
numbers, thus diffusion dominates. In the advection-dominated regime,
different approximations to $\widetilde{S}_{T}^{-1}$ are preferrable. In
contrast, CPR is designed for advection-dominated transport, and will likely
suffer for low Peclet numbers. The non-isothermal flow model problem is thus
an example where the optimal solution strategy depends on the solution regime.

##  5 Numerical experiments

To investigate the performance of the solver selection algorithm, we consider
two test cases: The Mandel problem in poromechanics, based on the model
presented in Section 4.1, and a case of non-isothermal flow and transport, as
described in Section 4.2. Full descriptions of the modeled systems are given
in A and B, respectively.

For these model problems we consider a series of experiments. In Section 5.1,
we consider cases with few candidate solver configurations to choose from, but
their optimal parameters must still be determined. In the poromechanical
simulation, we fix the categorical choices and focus only on optimizing the
numerical value of the stabilization parameter $L$ described in Section 4.1.
For the flow and transport problem, we seek when to switch between the CPR and
the Schur-complement based preconditioners from Section 4.2. This problem is
set up with a cyclic injection pattern, which repeatedly brings the problem
back to similar states, thus probing the selection algorithm’s ability to
exploit previous data. The cycle pattern is also constructed, by a somewhat
contrived design of the injection rate and schedule, to make the system
exhibit a wide range of Peclet numbers through the simulation, spanning
regimes where both the CPR and Schur-complement based preconditioner have
their forte. In Section 5.2, we consider the non-isothermal flow problem based
on the same setup as before but extend the number of candidate solver
configurations. This will test the selection algorithm’s ability to
efficiently zoom in on the best solvers.

As previously indicated, there are already well-established solver strategies
for both model problems, and we do not expect to arrive at results that are
significantly different from this existing knowledge. Indeed, our motivation
for considering these cases is, in part, to study how efficiently the solver
selection algorithm settles on strategies that could have been surmised from
existing understanding of these problems but without relying on heuristics and
manual tuning.

The solver selection framework is implemented in Python and is fully open-
source111github.com/Yuriyzabegaev/solver_selector. Moreover, we provide the
Docker image with the reproducible experiments from this article [48]. The
simulation models are implemented in the open-source package PorePy [49],
while the gradient boosting and Gaussian Process algorithms are provided by
Scikit-learn [50]. We use the implementation of GMRES and a direct solver
provided by SciPy [51], as well as the flexible GMRES (FGMRES) implementation
in PyAMG [52]. As preconditioners, we use the ILU(0) algorithm provided by
PETSc [53, 54] and the Boomer AMG library [55]. While other implementations
may be more efficient, our motivation is to select among the provided options;
thus, the important factor is to give the same settings to all options.

Linear systems are solved by GMRES with a restart after 30 inner iterations
and with a maximum of 10 restarts. The linear solver relative tolerance is
always set to $10^{-10}$. For the gradient boosting algorithm and Boomer AMG,
we use default parameters; parameters used for the Bayesian and heuristic
optimization methods are discussed in Section 3.

To limit measurement noise, we run our experiments with no other CPU- or
memory-intensive processes being executed on the same machine and use a single
core. We present the average and variance of 20 repetitions of all
experiments, with the solver selection potentially taking different random
choices in individual experiments. Data is not shared between these repeated
experiments.

###  5.1 Proving the concept of automated solver selection

In this set of experiments, we provide a small space of solver configurations
to the solver selector, and show that solver selection improves simulation
runtime compared to human-selected benchmarks. The optimization problem is
solved using the heuristic approach defined in Section 3.4 with the parameters
$\varepsilon=0.5$ and $\gamma=0.9$.

####  5.1.1 Optimization of numerical solver configuration parameter

To study selection of a numerical solver configuration parameter, we consider
Mandel’s problem and optimize the $L$ stabilization parameter of the Fixed-
Stress preconditioner. The solver space for this problem is:

$\mathcal{A}=\left\\{\text{GMRES}-\text{Fixed-
Stress}-L\in\left[L_{\text{min}},L_{\text{phys}}\right]-\text{Direct
subsolvers}\right\\}.$

The optimal value of $L$ for this specific simulated model is believed to be
static, so the context space $\mathcal{C}$ is void.

Figure 2(a) shows how the number of GMRES iterations depends on the value of
$L$ employed for a time step. We apply a linear scaling to $L$ so that
$L_{\text{min}}=0$ and $L_{\text{phys}}=1$, and the values from all the time
steps are shown together, so the results for the same $L$ can slightly vary
from one time step to another. From the figure, it is clear that the optimal
value for this problem is approximately 0.6. The search for the optimal value
of $L$ for one of the 20 runs of this experiment is shown in Figure 2(b): The
exploration process starts with random sampling and evolves to primarily
selecting the optimal value while reducing random decisions, following the
epsilon-greedy exploration policy. Continued exploration after the optimal
value has been found, and the associated overhead in computational cost, can
be seen as the price for safeguarding against dynamic changes in the optimal
value.

(a)

(b)

Figure 2: Optimization of numerical solver configuration parameter for the
Mandel’s problem. (a) Variation of number of GMRES iterations with the
stabilization coefficient $L$. (b) Evolution in $L$ from the solver selection
framework during one simulation.

To measure the overall computational gains from applying the solver selection
process, we consider the cumulative time spent on solving linear systems. In
Figure 3(a), this is shown for four solver options: The solver selection
algorithm, $L$ fixed to, respectively $L_{\text{min}}$ and $L_{\text{max}}$,
and $L$ randomly chosen in each time step in the interval
$[L_{\text{min}},L_{\text{max}}]$. As can be seen, the solver selection
delivers the best results, while the relative success of the randomly chosen
parameter can be contributed to the convex shape of the performance data shown
in Figure 2(a). The latter point indicates that the test problem is relatively
easy, but we reemphasize that the selection algorithm found the optimal $L$
value with no need for manual tuning, and with little overhead compared to
running the full simulation with the optimal value of $L$.

(a) Mandel’s problem

(b) Non-isothermal flow problem

Figure 3: Cumulative performance of solver configurations picked by the solver
selector compared to pre-selected configurations and random decisions. Lower
is better. The bold line of each color displays the mean value, while the
lighter-shaded region of the same color around it indicates the range from the
lowest to highest value among all the repeats of an experiment.

####  5.1.2 Optimization of categorical solver configuration choices

To study the impact of categorical decisions, we consider the non-isothermal
flow problem, with the following two solver configurations:

$\mathcal{A}=\left\\{\left(\text{GMRES}-\text{CPR}\right);\left(\text{GMRES}-\text{Schur}\right)\right\\}.$

All the numerical parameters are fixed. The Schur complement preconditioner in
this experiment applies the full factorization of the block inverse Eq. (21).
As noted in Section 4.2, the optimal solver configuration is expected to
depend on the dominant transport mechanism, measured by the Peclet number.
Accordingly, the Peclet number is included in the context space $\mathcal{C}$
together with other information listed in Table 1.

Parameter | Transformations
---|---
Time step value, s | $\log(\cdot)$
Peclet number | $\log\text{max}(\cdot)$
Peclet number | $\log\text{mean}(\cdot)$
Injection source rate, $\text{m}^{3}\text{s}^{-1}$ | $\log(\cdot)$
Production source rate, $\text{m}^{3}\text{s}^{-1}$ | $\log(\cdot)$
Well activity binary indicator |
Table 1: The non-isothermal flow problem context space. The right column
describes the transformations applied to these parameters.

The simulation goes through three injection periods that change the simulation
characteristics. We can track the process of solver selection in Figure 4,
which shows the dependence of the time step magnitude and the volumetric
injection rate on the time step number. The time step size increases rapidly
in periods without injection, followed by rapid decrease when injection is
switched on. If the Newton solver reaches the maximum number of iterations, we
decrease the time step and select a solver configuration again. The
performance data of these unsuccessful iterations is also used to update the
machine learning algorithm, and they are also included in our reporting on
linear solver behavior.

(a) Evolution of time step magnitude.

(b) Evolution of inlet source rate magnitude.

Figure 4: Solver selection for the non-isothermal flow problem. Different
colors represent the solver which was used for each time step.

Figures 4 further show the choice of solver configuration for each time step.
During the first injection period, the solver selector prioritizes
exploration, with frequent switching between the CPR and Schur complement
options. After approximately time step 30, the solver has collected sufficient
data to settle on the following policy: CPR is preferred when the injection
rate is high or the injection is switched off, while the Schur method is
chosen for low injection rates. Periods of high and low injection rates
correspond to high and low Peclet numbers, thus the transition from CPR to
Schur is as expected, and in agreement with the results from [16]. When
injection is switched off, the system is diffusion dominated, but it is also
close to steady state, thus the linear system is easier to solve. In this
regime, the ILU method applied in CPR is sufficiently efficient to result in a
lower total computational cost than that of the Schur method, although upon
inspection, the difference turns out to be minor.

To evaluate the effect of solver selection for this problem, we compare the
performance of the simulation driven by solver selection with three
benchmarks: Two preconditioners that use, respectively, the CPR and Schur
complement preconditioners throughout the simulation, and one which randomly
picks among them for each time step. The cumulative time to solve linear
systems is shown in Figure 3(b). For this specific problem setup, the relative
cost of misusing CPR is high, while the relative cost of misusing Schur is
relatively low. The extra cost of CPR is to a large degree associated with
three short time periods which coincide with the injection well switching from
rate to pressure control. The solver selector switches to the Schur complement
method at these points, thus avoiding the problem. When also accounting for
the overhead associated with exploration, the solver selector brings only
marginal gains compared to the Schur complement preconditioner. However,
solver performance will depend on the problem at hand, thus the optimal solver
cannot be known a priori. Moreover, the example also shows how the optimal
solver may change during a simulation, and this effect is expected to be
stronger if the problem complexity or range of solver options is increased.
The solver selector’s ability to adapt to such changes gives it the potential
to significantly reduce the computational burden of such simulations.

####  5.1.3 Exploration sensitivity

To better understand the behavior of the heuristic optimization algorithm, we
investigate its sensitivity to the exploration parameter $\varepsilon$. For
both the poroelastic and non-isothermal flow setups, we repeated the
experiments with a default value $\varepsilon=0.5$ (same as in the previous
experiment), with a higher value $\varepsilon=0.7$ and with zero exploration,
$\varepsilon=0$. The latter case corresponds to a simple greedy algorithm
which does not actively seek new data points, but may still obtain new data as
part of the optimization loop.

The results are shown in Table 2, where we also provide the results of the
random choice policy as a reference. For the poroelastic model, neglecting
exploration can significantly increase the spread of the results, as the
method may or may not approach the minimum of the convex function to be
optimized, cf. Figure 2(a). Increased exploration gives insignificant changes
to the standard model. The structure of the non-isothermal flow problem is
different, as its dynamic nature means the simulation characteristics change
for every time step. Hence, during the first injection period, the machine
learning algorithm always finds itself in an unexplored region of the space
$\mathcal{C}$. Predictions by machine learning algorithm outside the training
data region are unreliable and chaotic, so it naturally leads to a forced
exploration. Additional algorithmic exploration triggered by the $\varepsilon$
parameter prevents the worst-case result, and its overhead is negligible
compared to the mean- and best-case results of the greedy algorithm. However,
the amplified exploration leads to a worse result compared to the standard
algorithm because, for this simulated model, the cost of suboptimal decisions
is higher.

While these examples indicate that tuning of optimization parameters can give
some computational savings, they also show that the impact of such changes is
highly problem dependent. The main conclusion to be drawn from Table 2 is
that, while too much exploration in some cases can be considered an
unnecessary overhead, there may also be considerable costs in not exploring.

| Mandel’s problem | Non-isothermal flow problem
---|---|---
| Worst | Mean | Best | Worst | Mean | Best
Standard exploration $\varepsilon=0.5$ | 71.68 | 69.40 | 67.61 | 122.75 | 114.87 | 107.93
No exploration $\varepsilon=0$ | 102.42 | 70.19 | 65.65 | 134.64 | 114.57 | 107.83
More exploration $\varepsilon=0.7$ | 71.30 | 69.73 | 68.17 | 135.63 | 117.41 | 108.23
Random choice | 88.47 | 85.48 | 83.03 | 202.02 | 182.64 | 143.52
Table 2: Cumulative time (seconds) to solve linear systems for the whole
simulation for the Mandel’s problem and the non-isothermal flow problem.
Results for the solver selection policy with different values of the
$\varepsilon$ parameter and the result for the random policy. The minimal
value in each column is given in bold.

####  5.1.4 Knowledge transfer

The simulation context allows for the transfer of performance data between
simulations with different setups. For this to be meaningful, the relative
performance of different solvers must be similar in the setups, and this
similarity must be seen by the simulation context. Borrowing from machine
learning terminology, we refer to such knowledge transfer as warm-start
simulations, in contrast to the cold-start simulations where all performance
data must be collected through online learning.

To illustrate potential gains of warm starts, we alter the non-isothermal flow
setup by moving the two wells to regions of lower permeability, see Figure 8
in B. For both warm and cold start we consider 20 simulations to chart the
impact of random variations. The warm-start simulations are given access to
the data set $D$ from one of the experiments reported above. The results are
reported in Table 3, where we note that, since we are solving a problem with
different well locations, the time measurements cannot be directly compared to
those listed in Table 2. The warm start consistently leads to slightly better
results. While this may be considered a modest gain, we recall from Figure
4(a) that the initial phase of exploration, which is where the access to
previous data can make a difference compared to the cold start, only lasted
about 30 time steps with the original well locations and is likely of similar
length in this case. Again, this effect will grow with the problem complexity
or number of solver options.

| Worst | Mean | Best
---|---|---|---
Cold start | 59.17 | 55.42 | 50.14
Warm start | 56.99 | 54.19 | 50.06
Table 3: Cumulative time (seconds) to solve linear systems for the whole
simulation for the non-isothermal flow problem with the changed well
locations. A minimal value in each column is in bold. The time measurements
cannot be directly compared to those listed in Table 2 because we are solving
a problem with different well locations.

###  5.2 Extension of the available solver configuration space

Having studied the behavior of the solver selector applying the heuristic
optimization on a very limited set of choices, we now introduce two
extensions, focusing on the non-isothermal flow setup in Section 5.1.2. First,
we consider a wider set of categorical choices with in total 19 different sets
of options. This expansion will in effect lengthen the exploration phase, in
which the solver selector must deal with substantial uncertainty in the
simulation data. Our second extension is therefore to consider optimization
based on Gaussian Processes in addition to the heuristic approach, with the
motivation that GP’s explicit representation of uncertainty may make a
difference in the rate at which a policy is settled upon.

Id | Outer solver | Preconditioner | Factorization | Thermal | Flow
---|---|---|---|---|---
1 | GMRES | Schur | Lower | AMG | AMG
2 | GMRES | Schur | Upper | AMG | AMG
3 | GMRES | Schur | Full | AMG | AMG
4 | GMRES | CPR |  |  | AMG
5 | FGMRES | Schur | Lower | GMRES + AMG | GMRES + AMG
6 | FGMRES | Schur | Lower | AMG | GMRES + AMG
7 | FGMRES | Schur | Lower | GMRES + AMG | AMG
8 | FGMRES | Schur | Lower | AMG | AMG
9 | FGMRES | Schur | Upper | GMRES + AMG | GMRES + AMG
10 | FGMRES | Schur | Upper | AMG | GMRES + AMG
11 | FGMRES | Schur | Upper | GMRES + AMG | AMG
12 | FGMRES | Schur | Upper | AMG | AMG
13 | FGMRES | Schur | Full | GMRES + AMG | GMRES + AMG
14 | FGMRES | Schur | Full | AMG | GMRES + AMG
15 | FGMRES | Schur | Full | GMRES + AMG | AMG
16 | FGMRES | Schur | Full | AMG | AMG
17 | FGMRES | CPR |  |  | AMG
18 | FGMRES | CPR |  |  | GMRES + AMG
19 | Direct |  |  |  |
Table 4: The solver configurations considered in Section 5.2. Items that are
not relevant for a given configuration are blank.

The different solver configurations are listed in Table 4 and can be
summarized as follows: In addition to solving the diffusive subproblems with
AMG, as was done in Section 5.1.2, we also consider AMG preconditioned with
GMRES as inner solvers. This requires changing the outer Krylov solver to
flexible GMRES (FGMRES). For the preconditioners considered in Section 5.1.2,
i.e., without an inner GMRES method, we consider both GMRES and FGMRES as the
outer Krylov method. Recall that, while these methods are identical in terms
of the number of iterations, we consider implementations from two different
libraries. Thus, this is an example of how the solver selector can deal with
choices relating to software. For the Schur complement method, we consider the
three approximations of the inverse Jacobian Eq. (21) noted in Section 4.2.
Finally, we add the option of solving the entire system with a direct solver.

(a) Heuristic

(b) Gaussian Process

Figure 5: Histograms of the invocations of different solver configurations.
Average values and standard deviations are taken w.r.t. 20 repeats of the
simulation. The numbering of the options corresponds to that given in Table 4.
Figure 6: Time series of solver options and time step size for the heuristic
and GP optimization approaches, shown for a single simulation. The solvers are
grouped into five main categories defined in Section 5.2.

(a)

(b)

Figure 7: (a) Cumulative performance of solver configurations picked by the
solver selector in Experiment 5.2 for the heuristic and GP approaches (dashed
lines). The data is compared with the results of the Experiment 5.1.2 (solid
lines). The bold line of each color displays the mean value, while the
lighter-shaded region of the same color around it indicates the range from the
lowest to highest value among all the repeats of an experiment. (b) Overheads
to train the machine learning model for the heuristic and GP approaches, shown
for a single simulation and the instance of the machine learning model most
commonly invoked.

Figure 5 shows a histogram of the number of invocations of the different
solver configurations for the heuristic and GP optimization approaches. The
data is expressed in terms of the mean and standard deviation over 20 repeated
runs. A comparison of the results for the heuristic and GP optimization
methods indicates that the latter rapidly settles on a few preferred choices,
while the former distributes its choices more evenly among the available
options. The heuristic method has a clear preference for the Schur complement
approach, possibly reflecting that there are five times as many Schur options
than CPR, thus the former will be sampled more under random exploration. The
GP optimization leads to the Schur and CPR preconditioners being called
approximately equally often, and a comparison between solvers 4 and 17
(referring to Table 4) indicates a slight preference for the implementation
behind FGMRES over that of GMRES. Furthermore, the direct solver was hardly
chosen by any of the optimization methods.

To better understand the choices made, Figure 6 shows data from a single run
in the same experiment, with the options grouped into five clusters: CPR and
Schur with and without an inner GMRES method, and a direct solver. The GP
optimization seems to settle on a mixture of CPR and the Schur complement
approach reminiscent of the optimal strategy from Section 5.1.2, although the
policy is less clear with the higher number of options. In contrast, the
heuristic method hardly considers CPR-based options at all.

The result of these policies in terms of the time spent solving linear systems
is depicted in Figure 7(a). Despite having to experiment with 19 different
solvers, the average performance obtained by using GP optimization is
comparable to that of CPR, and in the best case, GP approaches the pure Schur
complement method. The heuristic optimization fails to find an efficient
policy and makes suboptimal choices for most time steps. A comparison with
Figure 3(b) shows that the variations among the 20 independent runs is
considerably larger than when choosing between only two options. This holds
for both the heuristic and GP optimization, and is a natural consequence of
the larger number of random decisions being made during exploration.

While the GP optimization produces superior results in terms of time spent on
linear solvers, this comes to the price of a significantly higher cost of the
optimization problem itself. Figure 7(b) depicts the time spent updating the
optimization model as new data is provided after each time step. For both the
heuristic and GP methods, we show measurements for the categorical choices
that were invoked most frequently, and note that these may be different
between the two optimizations. The cost of updating the heuristic model is
almost constant and low, while updating the GP model is significantly more
expensive, and grows with the size of the data to be fit. We attribute the
varying growth rate for the GP to the implementation of the training process,
which iterative searches for the correct algorithm parameters discussed in
Section 3.3. The data for the other algorithmic choices is similar, although
the growth in the cost of GP is less pronounced since these options have been
chosen less often. Thus, for this particular example, the total cost of
solving linear systems and maintaining the optimization method is lower for
the heuristic method than for the GP. However, this observation cannot be
readily transferred to other cases, as the cost of the two components scale
with different factors: The time spent on the linear solver depends on the
system size, coupling strength, etc., while the GP update depends on the
number of data points. The cost of the GP update can be partly alleviated, for
instance, by feeding it data less frequently and in batches, or by more
advanced techniques; see e.g., [56, 40]. Thus, the decision on which
optimization method to use depends on the simulation setup and technology
available.

##  6 Conclusion

We have presented a framework for selecting linear solvers in time-dependent
multiphysics simulations. The framework chooses from a set of user-provided
solver configurations that may include splitting schemes, preconditioners,
Krylov subspace methods and numerical parameters thereof. The selection
problem is formulated as a black-box optimization problem, which we solve by
Bayesian optimization with Gaussian Processes or a heuristic approach. The key
feature of our framework is that it dynamically adapts to performance data
generated during the simulation and may switch solver configuration between
time steps based on this information. This allows for quick determination of
promising solver configurations without relying on pre-generated training
data. Our framework may also incorporate characteristics of the problem and
the solution to dynamically adapt the linear solver.

The proposed framework was applied to two model problems, involving
poromechanics and non-isothermal single-phase fluid flow in porous media, for
both of which efficient linear solvers are known. Our framework could
successfully find efficient solver configurations for both problems and
perform better than the pre-selected benchmarks. The simulations also
illustrated how the solver selection problem becomes harder as the number of
candidate solvers grows and thus highlighted the utility in combining our
framework with existing knowledge of what are efficient solver configurations.
Finally, comparison between the two optimization methods showed that the
Gaussian Process gave the most accurate predictions of solver performance;
however, it also had a significantly higher cost than the heuristic approach,
thus the choice of optimization method is in itself non-trivial.

The task of selecting linear solvers for complex simulations is at the same
time critical for performance and difficult, in particular for users of
simulation software who have limited knowledge of such methods. As our
optimization approach automatizes this task, the cost and risk associated with
adding options to the pool of candidate solvers is reduced, thus enabling
experimentation with new and promising algorithms or implementations. Our
methodology can be inserted into existing simulation frameworks with minimal
adaptation and we therefore believe it can be a useful tool to speed up
numerical computations.

## Acknowledgment

This project has received funding from the VISTA program, The Norwegian
Academy of Science and Letters and Equinor and from the Norwegian Research
Council, Grant 308733.

##  Appendix A Model setup for poromechanics problem

The Mandel’s problem is an established verification setup for the
consolidation problem [57, 49]. We solve the equation system (4) in a
rectangular 2D domain. The domain is covered by two planes from the north and
south boundaries. Initially, a normal load is applied on these sides. The east
and west boundaries are stress-free and allow fluid to drain. The vertical
displacement at the top of the domain is time-dependent and given by the exact
solution. The gravity term is neglected. The symmetry of the problem allows us
to focus on the positive quarter of the domain, see [57].

We run the Mandel’s problem with the constant time step $\Delta t=10$ s during
100 time steps. The values of the material parameters for the solid defined in
Section 4.1 are the following: $\mu_{r}=2.475\cdot 10^{9}\text{ Pa}$,
$\lambda=1.65\cdot 10^{9}\text{ Pa}$, $m^{-1}=6.0606\cdot 10^{-11}\text{
Pa}^{-1}$, $\mathbf{K}=9.869\cdot\mathbf{I}\cdot 10^{-14}\text{ m}^{2}$
(homogeneous isotropic permeability), $b=1$. For the fluid, $\rho=10^{3}\text{
kg}\cdot\text{m}^{-3}$, $\mu=10^{-3}\text{ Pa}\cdot\text{s}$. The simulated
units are in SI; pressure is scaled to MPa.

##  Appendix B Model setup for non-isothermal flow in porous media

We solve the equation system (17) in a rectangular 2D domain of $(60\times
220)$ cells with heterogeneous porosity and permeability. The porosity and
permeability fields are taken from the first slice of the SPE10 benchmark
[58], model 2. To make permeability isotropic, we take only the first diagonal
component of the permeability tensor. The permeability field is multiplied by
1000 to achieve more hyperbolic properties of heat transport. The effect of
gravity is neglected. The model uses SI units, but the pressure is scaled to
MPa to improve the conditioning of the system.

The values of the material parameter for the solid are the following:
$\rho_{r}=2650\text{ kg}\cdot\text{m}^{-3}$, $c_{r}=920\text{
J}\cdot\text{kg}^{-1}\cdot\text{K}^{-1}$, $k_{Tr}=1.7295772056\text{
W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}$. For the fluid, $c_{f}=2093.4\text{
J}\cdot\text{kg}^{-1}\cdot\text{K}^{-1}$.

The initial pressure in the model is $41.369$ MPa and the initial temperature
is $288.706$ K. The fluid in the model is water. All boundary conditions
forbid fluid flow and energy transfer through the boundaries. The model
contains two wells: one for injection and one for production. A driving force
of the evolution in the model is a cyclic injection process. The model runs
for 300 simulated years. Three injection periods occur during the simulation:
at $t=0,50$ and $100$ years. Each injection period lasts 100 days; the hot
fluid is injected at a temperature of $422.039$ K, a maximum injection
pressure is $68.95$ MPa, and a maximum injection flux is $10^{-3}$ m3 s-1. The
production well is limited by a minimum pressure of $27.579$ MPa and a maximum
flux of $10^{-3}$ m3 s-1. The rest of the time, both wells are disabled.

Figure 8: Well locations and permeability of the non-isothermal fluid model in
4.2.

The setup was significantly motivated by the experiment from [16], where the
authors consider the Schur complement approximation preconditioner and the CPR
preconditioner, described in 4.2. One of the results shows that the Schur
complement approximation preconditioner performs better for the original SPE10
permeability, exhibiting more parabolic behavior. On the other hand, the CPR
preconditioner performs better for the permeability multiplied by 1000,
exhibiting more hyperbolic behavior. We model cyclic injection to witness both
solution regimes within one simulation to illustrate the case when dynamic
switching of solvers is beneficial. All the constitutive laws and the material
parameters are the same as in the numerical experiments in [16].

#### Sources

We consider the source/sink terms that represent injection and production
wells. The well locations are shown in Figure 8. The primarily used simulated
model utilizes well 0 for injection and well 1 for production. The alternative
simulated model in Section 5.1.4 utilizes well 1 for injection and well 2 for
production. A simple way to model these is by using point sources/sinks:

$f(\mathbf{x})=\sum_{i}q^{i}_{\text{inj}}(p,T)\delta(\mathbf{x-\mathbf{x^{i}_{\text{inj}}}})\rho(p,T_{\text{inj}})-\sum_{j}q^{j}_{\text{prod}}(p,T)\delta(\mathbf{x}-\mathbf{x}^{j}_{\text{prod}})\rho(p,T),$
(23)
$f_{T}(\mathbf{x})=\sum_{i}q^{i}_{\text{inj}}(p,T)\delta(\mathbf{x-\mathbf{x^{i}_{\text{inj}}}})\rho(p,T_{\text{inj}})cvT_{\text{inj}}-\sum_{j}q^{j}_{\text{prod}}(p,T)\delta(\mathbf{x}-\mathbf{x}^{j}_{\text{prod}})\rho(p,T)c_{v}T,$
(24)

where $\mathbf{x_{\text{inj}}}$ and $\mathbf{x}_{\text{prod}}$ represent the
location of injection and production wells, respectively. The Dirac delta
$\delta(\mathbf{x})$ is used to denote the source / sink locations,
$q^{i}_{\text{inj}}$ and $q_{\text{prod}}^{j}$ are injection and production
volumetric rates of the wells, respectively. We consider the Peaceman well
model [46, 47]. In our case of isotropic media, the rates are given by:

$q=\dfrac{2\pi h_{w}K}{\mu\ln(r_{e}/r_{w})}(p_{bh}-p),$ (25)

where $h_{w}$ is the height of well opening, $K$ is the isotropic scalar
permeability, $p_{bh}$ is the bottom-hole pressure, $r_{w}$ is the well
radius, and $r_{e}$ is the equivalent radius which is calculated by using the
formula:

$r_{e}=0.14\sqrt{D_{x}^{2}+D_{y}^{2}},$ (26)

where $D_{x}$ and $D_{y}$ are the horizontal lengths of the grid cell. We
choose $h_{w}=5$ meters and $r_{w}=0.1$ meters.

#### Density

The following formula is used for density:

$\rho(p,T)=\rho_{0}e^{c\left(p-p_{0}\right)}e^{\beta\left(T-T_{0}\right)},$
(27)

where $\rho_{0}$, $p_{0}$ and $T_{0}$ are reference values, $c$ is a
compressibility coefficient and $\beta$ is a thermal expansion coefficient.
The values of the parameters are representative of those used in reservoir
simulations: $p_{0}=1.01325$ bar, $T_{0}=228.7056$ K, $c=5.5\times 10^{-5}$
bar-1, $\beta=2.5\times 10^{-4}$ K-1 and $\rho_{0}=999$ kg m-3.

#### Viscosity

Viscosity is determined by the following correlation [59]:

$\mu(T_{F})=10^{A_{1}\gamma_{API}+A_{2}}T_{F}^{A_{3}\gamma_{API}+A_{4}},$ (28)

which takes temperature $T_{F}$ in ${}^{\circ}F$ and returns viscosity in cp
(0.001 kg m-1 s-1). The values of dimensionless parameters are given:
$A_{1}=-0.8021$, $A_{2}=23.8765$, $A_{3}=0.31458$, $A_{4}=-9.21592$. The value
of $\gamma_{API}$ is 10 (for water).

## References

  * Settari and Mounts [1998] A. Settari and F.M. Mounts.  A Coupled Reservoir and Geomechanical Simulation System.  _SPE J._ , 3(03):219–226, 09 1998.  ISSN 1086-055X.  doi:10.2118/50939-PA.
  * White et al. [2016] J.A. White, N. Castelletto, and H.A. Tchelepi.  Block-partitioned solvers for coupled poromechanics: A unified framework.  _Comput. Methods Appl. Mech. Engrg._ , 303:55–74, 5 2016.  ISSN 00457825.  doi:10.1016/j.cma.2016.01.008.
  * Both et al. [2017] J.W. Both, M. Borregales, J.M. Nordbotten, K. Kumar, and F.A Radu.  Robust fixed stress splitting for biot’s equations in heterogeneous media.  _Appl. Math. Lett._ , 68:101–108, 2017.  doi:10.1016/j.aml.2016.12.019.
  * Kim et al. [2011a] J. Kim, H.A. Tchelepi, and R. Juanes.  Stability and convergence of sequential methods for coupled flow and geomechanics: Drained and undrained splits.  _Comput. Methods Appl. Mech. Engrg._ , 23(200):2094–2116, 2011a.  doi:10.1016/j.cma.2011.02.011.
  * Kim et al. [2011b] J. Kim, H.A. Tchelepi, and R. Juanes.  Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits.  _Comput. Methods Appl. Mech. Engrg._ , 13(200):1591–1606, 2011b.  doi:10.1016/j.cma.2010.12.022.
  * Mikelić and Wheeler [2013] A. Mikelić and M.F. Wheeler.  Convergence of iterative coupling for coupled flow and geomechanics.  _Comput. Geosci._ , 17:455–461, 6 2013.  ISSN 1420-0597.  doi:10.1007/s10596-012-9318-y.
  * Storvik et al. [2018] E. Storvik, J.W. Both, K. Kumar, J.M. Nordbotten, and F.A. Radu.  On the optimization of the fixed-stress splitting for biot’s equations.  _Int. J. Numer. Methods Eng._ , 120:179–194, 11 2018.  doi:10.1002/nme.6130.
  * Trangenstein and Bell [1989] J.A. Trangenstein and J.B. Bell.  Mathematical structure of compositional reservoir simulation.  _SIAM J. Sci. Stat. Comput._ , 10(5):817–845, 1989.  doi:10.1137/0910049.
  * Wallis [1983] J.R. Wallis.  Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration.  In _Proc. of SPE Reservoir Simulation Symposium_ , pages SPE–12265–MS, 1983.  doi:10.2118/12265-MS.
  * Wallis et al. [1985] J.R. Wallis, R.P. Kendall, and T.E. Little.  Constrained Residual Acceleration of Conjugate Residual Methods.  In _Proc. of SPE Reservoir Simulation Symposium_ , pages SPE–13536–MS, 1985.  doi:10.2118/13536-MS.
  * Lacroix et al. [2003] S. Lacroix, Yu. Vassilevski, J. Wheeler, and M. Wheeler.  Iterative solution methods for modeling multiphase flow in porous media fully implicitly.  _SIAM J. Sci. Comput._ , 25(3):905–926, 2003.  doi:10.1137/S106482750240443X.
  * Stüben et al. [2007] K. Stüben, T. Clees, H. Klie, Bo Lu, and M.F. Wheeler.  Algebraic Multigrid Methods (AMG) for the Efficient Solution of Fully Implicit Formulations in Reservoir Simulation.  In _Proc. of SPE Reservoir Simulation Symposium_ , pages SPE–105832–MS, 2007.  doi:10.2118/105832-MS.
  * Cao et al. [2005] H. Cao, H.A. Tchelepi, J. Wallis, and H. Yardumian.  Parallel Scalable Unstructured CPR-Type Linear Solver for Reservoir Simulation.  In _SPE Annual Technical Conference and Exhibition_ , volume All Days, pages SPE–96809–MS, 10 2005.  doi:10.2118/96809-MS.
  * Gries et al. [2014] S. Gries, K. Stüben, G.L. Brown, D. Chen, and D.A. Collins.  Preconditioning for efficiently applying algebraic multigrid in fully implicit reservoir simulations.  _SPE J._ , 19(04):726–736, 2014.  doi:10.2118/163608-PA.
  * Li et al. [2015] G. Li, J. Wallis, and G. Shaw.  A Parallel Linear Solver Algorithm for Solving Difficult Large Scale Thermal Models.  In _Proc. of SPE Reservoir Simulation Symposium_ , pages SPE–173207–MS, 2015.
  * Roy et al. [2019] T. Roy, T.B. Jönsthövel, C. Lemon, and A.J. Wathen.  A block preconditioner for non-isothermal flow in porous media.  _J. Comput. Phys._ , 395:636–652, 10 2019.  ISSN 00219991.  doi:10.1016/j.jcp.2019.06.038.
  * Cremon et al. [2020] M.A. Cremon, N. Castelletto, and J.A. White.  Multi-stage preconditioners for thermal–compositional–reactive flow in porous media.  _J. Comput. Phys._ , 418:109607, 2020.  doi:10.1016/j.jcp.2020.109607.
  * Roy et al. [2020] T. Roy, T.B. Jönsthövel, C. Lemon, and A.J. Wathen.  A constrained pressure-temperature residual (cptr) method for non-isothermal multiphase flow in porous media.  _SIAM J. Sci. Comput._ , 42(4):B1014–B1040, 2020.  doi:10.1137/19M1292023.
  * Rice [1976] J.R. Rice.  The Algorithm Selection Problem.  In _Advances in Computers_ , volume 15, pages 65–118. Elsevier, 1976.  ISBN 978-0-12-012115-1.  doi:10.1016/S0065-2458(08)60520-3.
  * K.A. [2009] Smith-Miles K.A.  Cross-disciplinary perspectives on meta-learning for algorithm selection.  _ACM Comput. Surv._ , 41:1–25, 1 2009.  ISSN 0360-0300.  doi:10.1145/1456650.1456656.
  * Bischl et al. [2016] B. Bischl, P. Kerschke, L. Kotthoff, M. Lindauer, Y. Malitsky, A. Fréchette, H. Hoos, F. Hutter, K. Leyton-Brown, K. Tierney, and J. Vanschoren.  Aslib: A benchmark library for algorithm selection.  _Artif. Intell._ , 237:41–58, 2016.  ISSN 0004-3702.  doi:10.1016/j.artint.2016.04.003.
  * Jessup et al. [2016] E. Jessup, P. Motter, B. Norris, and K. Sood.  Performance-based numerical solver selection in the lighthouse framework.  _SIAM J. Sci. Comput._ , 38(5):S750–S771, 2016.  doi:10.1137/15M1028406.
  * Weerawarana et al. [1996] S. Weerawarana, E.N. Houstis, J.R. Rice, A. Joshi, and C.E. Houstis.  PYTHIA: A knowledge-based system to select scientific algorithms.  _ACM Trans. Math. Softw._ , 22(4):447–468, December 1996.  ISSN 0098-3500.  doi:10.1145/235815.235820.
  * Dongarra et al. [2006] J. Dongarra, G. Bosilca, Z. Chen, V. Eijkhout, G. E. Fagg, E. Fuentes, J. Langou, P. Luszczek, J. Pjesivac-Grbovic, K. Seymour, H. You, and S. S. Vadhiyar.  Self-adapting numerical software (sans) effort.  _IBM J. Res. Dev._ , 50(2.3):223–238, March 2006.  ISSN 0018-8646.  doi:10.1147/rd.502.0223.
  * Bhowmick et al. [2011] S. Bhowmick, V. Eijkhout, Y. Freund, E. Fuentes, and D. Keyes.  Application of Alternating Decision Trees in Selecting Sparse Linear Solvers.  In Ken Naono, Keita Teranishi, John Cavazos, and Reiji Suda, editors, _Software Automatic Tuning_ , pages 153–173. Springer New York, New York, NY, 2011.  ISBN 978-1-4419-6934-7 978-1-4419-6935-4.  doi:10.1007/978-1-4419-6935-4_10.
  * Eller et al. [2012] P.R. Eller, J.C. Cheng, and R.S. Maier.  Dynamic linear solver selection for transient simulations using machine learning on distributed systems.  In _2012 IEEE 26th International Parallel and Distributed Processing Symposium Workshops & PhD Forum_, pages 1915–1924, May 2012.  doi:10.1109/IPDPSW.2012.239.
  * Liu et al. [2021] Y. Liu, W.M. Sid-Lakhdar, O. Marques, X. Zhu, C. Meng, J.W. Demmel, and X.S. Li.  Gptune: Multitask learning for autotuning exascale applications.  In _Proceedings of the 26th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming_ , PPoPP ’21, page 234–246. Association for Computing Machinery, 2021.  ISBN 9781450382946.  doi:10.1145/3437801.3441621.
  * Roy et al. [2021] R.B. Roy, T. Patel, V. Gadepally, and D. Tiwari.  Bliss: auto-tuning complex applications using a pool of diverse lightweight learning models.  In _Proceedings of the 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation_ , pages 1280–1295, Virtual Canada, June 2021. ACM.  ISBN 978-1-4503-8391-2.  doi:10.1145/3453483.3454109.
  * Rasmussen et al. [2006] C.E. Rasmussen, C. Williams, et al.  _Gaussian processes for machine learning_ , volume 1.  Springer, 2006.
  * Mishev et al. [2009] I. Mishev, B. Beckner, N. Fedorova, and S. Terekhov.  Linear Solver Performance Optimization in Reservoir Simulation Studies.  In _Proc. of SPE Reservoir Simulation Symposium_ , pages SPE–119266–MS, 2009.  doi:10.2118/119266-MS.
  * Bhowmick et al. [2006] S. Bhowmick, D. Kaushik, L. McInnes, B. Norris, and P. Raghavan.  Parallel adaptive solvers in compressible petsc-fun3d simulations*.  In A. Deane, A. Ecer, J. McDonough, N. Satofuka, G. Brenner, D.R. Emerson, J. Periaux, and D. Tromeur-Dervout, editors, _Parallel Computational Fluid Dynamics 2005_ , pages 277–284. Elsevier, Amsterdam, 2006.  ISBN 978-0-444-52206-1.  doi:10.1016/B978-044452206-1/50033-1.
  * Clees and Ganzer [2010] T. Clees and L. Ganzer.  An efficient algebraic multigrid solver strategy for adaptive implicit methods in oilreservoir simulation.  _SPE J._ , 15:670–681, 09 2010.  doi:10.2118/105789-PA.
  * Silva et al. [2021] V. L.S. Silva, P. Salinas, M. D. Jackson, and C. C. Pain.  Machine learning acceleration for nonlinear solvers applied to multiphase porous media flow.  _Comput. Methods Appl. Mech. Engrg._ , 384, 10 2021.  ISSN 00457825.  doi:10.1016/j.cma.2021.113989.
  * Jiranek et al. [2010] P. Jiranek, Z. Strakos, and M. Vohralík.  A posteriori error estimates including algebraic error: computable upper bounds and stopping criteria for iterative solvers.  _SIAM J. Sci. Comput._ , 32(3):1567–1590, 2010.  doi:10.1137/08073706X.
  * Ahmed et al. [2020] E. Ahmed, J. Nordbotten, and Florin Adrian Radu.  Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems.  _J. Comput. Appl. Math_ , 364:112312, 2020.  doi:10.1016/j.cam.2019.06.028.
  * Caldana et al. [2023] M. Caldana, P. F. Antonietti, and L.’ Dede.  A deep learning algorithm to accelerate algebraic multigrid methods in finite element solvers of 3d elliptic pdes, 2023.
  * Ramesh Kumar and Tene [2024] K. Ramesh Kumar and M. Tene.  Algebraic multiscale grid coarsening using unsupervised machine learning for subsurface flow simulation.  _J. Comput. Phys._ , 496:112570, 2024.  ISSN 0021-9991.  doi:10.1016/j.jcp.2023.112570.
  * Huang et al. [2023] R. Huang, R. Li, and Y. Xi.  Learning optimal multigrid smoothers via neural networks.  _SIAM J. Sci. Comput._ , 45(3):S199–S225, 2023.  doi:10.1137/21M1430030.
  * Nguyen et al. [2020] D. Nguyen, S. Gupta, S. Rana, A. Shilton, and S. Venkatesh.  Bayesian optimization for categorical and category-specific continuous inputs.  _Proc. Int. AAAI Conf._ , 34(04):5256–5263, Apr. 2020.  doi:10.1609/aaai.v34i04.5971.
  * Shahriari et al. [2016] B. Shahriari, K. Swersky, Z. Wang, R.P. Adams, and N. de Freitas.  Taking the human out of the loop: A review of bayesian optimization.  _Proc. IEEE_ , 104(1):148–175, 2016.  doi:10.1109/JPROC.2015.2494218.
  * Friedman [2001] J.H. Friedman.  Greedy function approximation: A gradient boosting machine.  _Ann. Stat._ , 29(5):1189 – 1232, 2001.  doi:10.1214/aos/1013203451.
  * Coussy [2004] O. Coussy.  _Poromechanics_.  John Wiley & Sons, 2004.
  * Aziz and Settari [1979] K. Aziz and A. Settari.  _Petroleum Reservoir Simulation_.  Springer Netherlands, 1979.  ISBN 9780853347873.
  * Nordbotten [2016] J.M. Nordbotten.  Stable cell-centered finite volume discretization for biot equations.  _SIAM J. Numer. Anal._ , 54(2):942–968, 2016.  doi:10.1137/15M1014280.
  * Both and Köcher [2019] J.W. Both and U. Köcher.  Numerical investigation on the fixed-stress splitting scheme for biot’s equations: Optimality of the tuning parameter.  In Florin Adrian Radu, Kundan Kumar, Inga Berre, Jan Martin Nordbotten, and Iuliu Sorin Pop, editors, _Numerical Mathematics and Advanced Applications ENUMATH 2017_ , pages 789–797, Cham, 2019. Springer International Publishing.  ISBN 978-3-319-96415-7.
  * Peaceman [1978] D.W. Peaceman.  Interpretation of Well-Block Pressures in Numerical Reservoir Simulation(includes associated paper 6988 ).  _SPE J._ , 18(03):183–194, 06 1978.
  * Chen and Zhang [2009] Z. Chen and Y. Zhang.  Well flow models for various numerical methods.  _Int. J. Numer. Anal. Mode._ , 6, 01 2009.
  * Zabegaev et al. [2024] Y. Zabegaev, E. Keilegavlen, E. Iversen, and I. Berre.  Runscripts for the manuscript ‘Automated solver selection for simulation of multiphysics processes in porous media’, January 2024.  doi:10.5281/zenodo.10407058.
  * Keilegavlen et al. [2021] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, and I. Berre.  Porepy: an open-source software for simulation of multiphysics processes in fractured porous media.  _Comput. Geosci._ , 25:243–265, 2 2021.  ISSN 15731499.  doi:10.1007/s10596-020-10002-5.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _J. Mach. Learn. Res._ , 12:2825–2830, 2011.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, İ. Polat, Y. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nat. Methods._ , 17:261–272, 2020.  doi:10.1038/s41592-019-0686-2.
  * Bell et al. [2023] N. Bell, L.N. Olson, J. Schroder, and B. Southworth.  PyAMG: Algebraic multigrid solvers in Python.  _J. Open Source Softw._ , 8(87):5495, 2023.  doi:10.21105/joss.05495.
  * Balay et al. [2023] S. Balay, S. Abhyankar, M.F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W.D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D Kaushik, M.G. Knepley, F. Kong, S. Kruger, May D.A., L. Curfman McInnes, Tran Mills R., L. Mitchell, T. Munson, J.E. Roman, K. Rupp, P. Sanan, J. Sarich, B.F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang.  PETSc/TAO users manual.  Technical Report ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023.
  * Dalcin et al. [2011] L.D. Dalcin, R.R Paz, P.A. Kler, and A. Cosimo.  Parallel distributed computing using Python.  _Adv. Water Resour._ , 34(9):1124 – 1139, 2011.  ISSN 0309-1708.  doi:10.1016/j.advwatres.2011.04.013.  New Computational Methods and Software Tools.
  * Henson and Meier Yang [2002] V.E. Henson and U. Meier Yang.  Boomeramg: A parallel algebraic multigrid solver and preconditioner.  _Appl. Numer. Math._ , 41(1):155–177, 2002.  ISSN 0168-9274.  doi:10.1016/S0168-9274(01)00115-5.  Developments and Trends in Iterative Methods for Large Systems of Equations - in memorium Rudiger Weiss.
  * Vakili et al. [2021] S. Vakili, H. Moss, A. Artemev, V. Dutordoir, and V. Picheny.  Scalable thompson sampling using sparse gaussian process models.  In M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , volume 34, pages 5631–5643. Curran Associates, Inc., 2021.
  * Mandel [1953] J. Mandel.  Consolidation des sols (Étude mathématique).  _Géotechnique_ , 3(7):287–299, 1953.  doi:10.1680/geot.1953.3.7.287.
  * Christie and Blunt [2001] M.A. Christie and M.J. Blunt.  Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques.  _SPE Reserv. Evaluation Eng._ , 4(04):308–317, 08 2001.  ISSN 1094-6470.  doi:10.2118/72469-PA.
  * Bennison [1998] T. Bennison.  Prediction of heavy oil viscosity.  In _IBC Heavy oil field development conference_ , volume 2, page 4. Citeseer, 1998.