In recent years, attention has grown around novel technologies and applications in the
subsurface, like geothermal energy production [1]}, [2]}, [3]}, hydraulic fracturing [4]}, [5]}, [6]}, CO\(_2\)  sequestration [7]}, [8]}, [9]} and underground gas storage [10]}, [11]}, [12]}. In these contexts, one of the key components is
the simultaneous simulation of frictional contact mechanics and fluid flow in faults and fractures, which represent tightly coupled
physical processes.
In fact, the aperture and slippage between the contact surfaces drive the fluid flow in the fractures, while the pressure variation perturbs the stress state in the surrounding medium and influences the contact mechanics itself.
To achieve the
desired accuracy, large domains are usually required, with high resolution representations
of geological structures and their heterogeneous properties [13]}, [14]}, and, specifically, of faults and
fracture networks [15]}, [16]}, [17]}, [18]}, [19]}, [20]}, [21]}, [22]}, [23]}, [24]}, [25]}, [26]}, [27]}, [28]}, [29]}, [30]}, [31]}, [32]}. It is, therefore, natural to have a growing demand towards the development of sophisticated models of increasing size, which are computationally intensive and require better and better performances. A key factor in this sense is the linear solver, which is usually by far the most
time-consuming component in a real-world simulation [33]}, [34]}.
