  Many subsurface engineering applications involve tight-coupling between fluid
flow, solid deformation, fracturing, and similar processes. To better
understand the complex interplay of different governing equations, and
therefore design efficient and safe operations, numerical simulations are
widely used. Given the relatively long time-scales of interest, fully-implicit
time-stepping schemes are often necessary to avoid time-step stability
restrictions. A major computational bottleneck for these methods, however, is
the linear solver. These systems are extremely large and ill-conditioned.
Because of the wide range of processes and couplings that may be involved--e.g.
formation and propagation of fractures, deformation of the solid porous medium,
viscous flow of one or more fluids in the pores and fractures, complicated well
sources and sinks, etc.--it is difficult to develop general-purpose but
scalable linear solver frameworks. This challenge is further aggravated by the
range of different discretization schemes that may be adopted, which have a
direct impact on the linear system structure. To address this obstacle, we
describe a flexible framework based on multigrid reduction that can produce
purely algebraic preconditioners for a wide spectrum of relevant physics and
discretizations. We demonstrate its broad applicability by constructing
scalable preconditioners for several problems, notably: a hybrid discretization
of single-phase flow, compositional multiphase flow with complex wells, and
hydraulic fracturing simulations. Extension to other systems can be handled
quite naturally. We demonstrate the efficiency and scalability of the resulting
solvers through numerical examples of difficult, field-scale problems.
