This paper deals with the numerical simulation of flow and transport in a fractured porous medium. We are particularly interested in the simulation of very large fracture networks (several thousands of fractures in a porous medium), and how numerical methods must be adapted to scale-up to such large systems.

We consider one phase Darcy flow of water, where ractures are represented as co-dimension one interfaces, and the conservation law in the fracture is modified to include a source term representing the exchange with the porous matrix. We assume the pressure is continuous across the matrix--fracture interface, and the model is closed by adding a zero total flux along fracture intersections.

A first step is to obtain realistic large scale networks. We rely on networks obtained from the DFN.Lab software (from teh Fractory) . A dedicated mesh generation code MODFRAC (a proprietary software owned by Inria and the University of Technology of Troyes) is then used to mesh the resulting network, and this is used as input to a 3D mesh generation software. 

The discretization is based on mixed-hybrid finite element. This choice presents several advantages: the methods is conservative by construction, and is robust with respect to heterogeneous or anisotropic media. Hybridization leads to a symmetric and positive definite matrix with one unknown on each face of the mesh , but its main interest is to give a very natural way of coupling the 3D and 2D models: the pressure at the center of triangular fracture elements coincides with the Lagrange multiplier (or pressure trace) on the faces or tetrahedra in the porous matrix.

The last, but crucial, step is to solve the resulting linear system. It is both large, because of the large number of unknowns in the porous matrix, and may be highly ill-conditioned because the 3D mesh is constrained by the geometry of the fractures, and may become quite distorted. Direct methods worked well for fracture networks (whose topology remains two-dimensional) or smaller examples but become impractical for large 3D meshes. One must turn to iterative methods, such as conjugate gradient or GMRES. The critical ingredient in ensuring the success of an iterative method is the choice of a preconditioner. We compare the performance of a sparse direct solver (16 hours CPU, 2TB memory) with that of a conjugate gradient with an algebraic multigrid preconditioner 7h CPU for 5736 iterations) , 8 GB memory). The iterative method reduces the memory usage, but the number of iterations remains very high.

An interesting alternative to algebraic multigrid is afforded by domain decomposition methods, in particular the recently developed multilevel spectral domain decomposition methods based on the GenEO framework. We report a preliminary performance result: on a model with 87 000 fractures, 16 million cells and 37 millions unknowns, AMG did not converge to a tolerance of 10^{-12} after 50 000 iterations, while the HPDDM preconditioner used from PETSc converged in 131 iterations, and 8 minutes on 1024 cores.