It is well known that iterative methods based on projections/orthogonalizations onto Krylov subspaces [1]} are in practice mandatory to solve large and sparse linear systems deriving from the discretization of
PDEs, because they allow for a lower complexity, smaller memory requirement, and better degree of
algorithmic parallelism than direct methods [2]}.
However, robustness, scalability and computational efficiency of this class of methods is tightly connected with the choice of a proper preconditioning technique [1]}. Roughly speaking, preconditioners are approximate applications of the system matrix inverse, and, from the algebraic viewpoint, can be
classified into three main categories: (i) incomplete factorizations [4]}, [5]}, [6]}, (ii) approximate inverses
[7]}, [8]}, [9]}, [10]}, [11]}, [12]}, and (iii)
multilevel methods, i.e., domain decomposition [13]}, [14]}, [15]}, [16]}, [17]} and multigrid-like techniques
[18]}, [19]}, [20]}, [21]}, [22]}, [23]}, [24]}, [25]}, [26]}, [27]}, [28]}, [29]}, [30]}, [31]}. A key feature for a modern preconditioning framework is the
algorithmic scalability, i.e., the ability to solve an increasingly refined problem with an approximately
constant number of iterations of the Krylov solver. This property is particularly important in view of the development of problems of increasing size by exploiting the availability of massively parallel computational platforms. Incomplete factorizations and approximate inverses can exhibit amazing performances, but do
not have a linear complexity with the system size. By distinction, multilevel methods can have a lower performance on a single system, but are designed to be optimal with respect to the scalability issue. Algebraic multigrid (AMG, [32]}) is one of the most effective multilevel
approaches and consists of the complementary use of: (i) a smoother that reduces high
frequency errors, (ii) a coarse grid correction that reduces low frequency errors, and (iii)
restriction and interpolation operators, to move from one grid to another. Starting from the original
works, e.g., [33]}, a wide range of multigrid approaches has appeared in the literature,
extending the applicability of this method, originally designed for elliptic PDEs, to both
non-symmetric [34]}, [35]} and block matrices [36]}, [37]}, [38]}, [39]}, [40]}, [41]}.
Nonetheless, robustness and efficiency is still an open issue for AMG whenever used as a black-box tool in problems with these algebraic properties.
The Jacobian matrix arising from the model considered herein is a
non-symmetric \(3 \times 3\)  block matrix and, despite the available studies for
similar problems, none of them can be straightforwardly and effectively
applied to our case. In the context of geomechanical simulations, only a few studies
on \(2 \times 2\)  block Jacobian systems [42]}, [43]}, [40]} are found by the authors.
