A hybrid high-order method for flow simulations in discrete fracture networks.

We are interested in solving ﬂow in large tridimensional Discrete Fracture Networks (DFN) with the hybrid high-order (HHO) method. The objectives of this paper are: (1) to demonstrate the beneﬁt of using a high-order method for computing macroscopic quantities, like the equivalent permeability of fracture rocks; (2) to present the computational eﬃciency of our C++ software, NEF++ , which implements the solving of ﬂow in fractures based on the HHO method


The flow problem in fractured rocks
In fractured rocks, fluid flows mostly within a complex arrangement of fractures, classically modeled as a Discrete Fracture Network (DFN) [1,2].In the present reduced model, the fractures, denoted by Ω f , f = 1, ..., N f , are distributed in a three-dimensional domain Ω and are modeled as ellipses whose position and orientation are evaluated from statistical laws given by geological studies [3,4].We consider single phase flow problems within these networks of fractures.As we are mainly interested in flow simulations in granite type rocks, a classical assumption is to consider the rock matrix as impervious.Figure 1 presents three examples of DFN in a cubic domain.Let x be the local 2D coordinates of fracture Ω f .Let N be the total number of intersections between fractures, I k be the k th intersection, k = 1, ..., N, and F k be the set of fractures containing I k .In each fracture Ω f , we assume that the governing equations for the hydraulic head scalar function p and for the flux per unit length function u are the mass conservation equation and Poiseuille's law [1]: The parameter T (x) is a given transmissivity field (unit [m 2 .s−1 ]).The function f ∈ L 2 (Ω f ) represents the sources/sinks.Additionally, continuity of the hydraulic head and continuity of the transversal flux apply at the intersections between the fractures [2,5]: where p k,i is the trace of hydraulic head on I k in fracture Ω f , p k is the unknown hydraulic head on the intersection I k and u k,i • n k, f is the normal flux through I k coming from fracture Ω f , with n k, f the outward normal unit vector of the intersection I k with respect to the fracture Ω f .Boundary conditions (BC) on the cube faces are of Dirichlet or Neumann type.For edges that belong to the border of the fractures but not to a cube face, a homogeneous Neumann BC is applied to express the imperviousness of the rock matrix.

The HHO method for solving flow in DFN
Several methods have been developed to solve flow in DFN in the recent years as detailed in the survey [6] and the references therein.The methods highly depend on the mesh strategy chosen to mesh the DFN [7,8].In our work, we keep the intersections explicitely.It implies a substantial work regarding the development of robust and efficient software in order to be able to mesh efficiently large networks with a good quality mesh [9].In all our test cases, the mesh is generated with the software BLSURF_FRAC [10,11] and the planar mesher is BL2D [12].The data files generated by BLSURF_FRAC follow the description given by Appendix A in [13].
Here we consider the so-called conforming discretizations at the intersections between the ellipses but the software BLSURF_FRAC is also able to generate non-conforming discretizations as well.The advantage of keeping the intersections explicitely in the mesh generation is that it allows to attach unknowns to the edges and then continuity conditions (2) are easier to impose.
Among the methods that attach unknowns to the edges, let us cite the mixed-hybrid finite elements method (MHFEM), for conforming [14,15,16,17,5] or non-conforming [18,19] discretizations at the intersections.More recently, a hybrid high-order (HHO) method has been developed [20,21].HHO is already used in many applications and has been recently used for fracture/matrix coupling [22].HHO is closely related to Hybridizable Discontinuous Galerkin methods (HDG) [23].The main advantages of HHO are (1) it allows general meshes (including polytopal cells and nonmatching interfaces), ( 2) it manages arbitrary polynomial face orders k, (3) it leads solve a linear system with only the unknowns at the edges and the matrix of this system is symmetric positive definite, (4) it delivers approximate solutions converging at order h k+1 in the energy norm and h k+2 in the L 2 -norm (if full elliptic regularity holds) [20,21].Moreover, this HHO method is implemented in the open source library, DiSk++ [24], which is a C++ template based library, both in the dimension and also in the finite element shapes.Notice that only the 2D feature of the DiSk++ library is used in this study as the rock matrix is assumed impervious, however the dimensional templating offered by DiSk++ will be very useful for future porous fractured rocks simulations.

Computation of the equivalent permeability with the HHO method
The goal of this section is to demonstrate the benefit of using a high-order method for computing upscaled quantities, like the equivalent permeability.
The equivalent permeability tensor is a macroscopic quantity of interest classically used by hydrogeologists for upscaling [25].Its components can be derived from numerical simulations.Typically, the three diagonal components of the permeability tensor are given by applying permeameter boundary conditions in the directions x, y and z respectively.As we are rather interested in analyzing the performance of the HHO method to compute such macroscopic quantities, we focus here only on a flow in the direction x to derive the x-component of the permeability tensor defined as: L ∆h for a cubic domain of size L, with Q in, x the input flux (units m 3 .s−1 ) with respect to permeameter boundary conditions in the direction x.
We propose to compute the equivalent permeability in the direction x of the small network B0 shown on Figure 2 (left).The domain is a cube of size L = 20.This network has 1, 397 fractures and 2, 481 intersections.We imposed a permeameter boundary condition in the direction x with a difference of hydraulic head of 10 m between the two opposite cube faces.The mean hydraulic head solution obtained with the MHFEM (Raviart-Thomas 0) for a mesh with 8, 533, 221 edges is shown on Figure 2 (right).We compute K x with the MHFEM RT0 3  The simulations for k = 1 and k = 2 on the finer mesh are not available yet as they require a lot of computational resources.For the MHFEM RT0 method, the software that is used is a Matlab software, called NEF-Flow [11], developed at Inria and CNRS (France), and for the HHO method for DFN, the simulations are performed with the C++ software NEF++, developed at Inria and described in more details in the following section.
Figure 3 presents the equivalent permeability with respect to the number of degrees of freedom (dofs).For the MHFEM RT0 and the HHO method with k = 0, the number of dofs are equal to the total number of edges, denoted by n E .For the HHO method with k = 1, the number of dofs are 2 n E .For the HHO method with k = 2, the number of dofs are 3 n E .At low order (k = 0), the curves for K x obtained with the HHO method and the MHFEM RT0 method almost superimpose, which was expected as the two methods are very close.The benefits of the increased orders of convergence are clearly seen from Figure 3 since the equivalent permeability is better approximated by a fixed number of dofs if resorting to a higher-order method.This result shows that the exact solution has enough regularity to take advantage of the computational efficiency delivered by higher-order methods.

Performance obtained with the NEF++ software
Solving the flow problem (1) − ( 2) within large 3D DFNs requires robust and efficient software.The goal of this section is to present the C++ software we have developed at Inria, called NEF++, and based on the C++17 standard.NEF++ relies on the Eigen library, which is a C++ template library for linear algebra [26] and on the DiSk++ library for HHO.The linear systems can be solved with direct solvers or with iterative solvers like the preconditioned conjugate gradient or multigrid solvers.In NEF++, the following two direct solvers can be called: Pardiso from Intel MKL library [27] or SuiteSparse [28].Both solvers support tasks parallelism, either using OpenMP or Intel TBB and support SIMD (Single Instruction, Multiple Data) vectorization.
We propose to solve flow in the three DFN B1, B2 and B3 shown on Figure 1.We imposed a permeameter boundary condition in the direction x with a difference of hydraulic head of 10 m between the two opposite cube faces.The transmissivity is taken as a constant per fracture and is different from one plane to another.We consider a cubic domain Ω of size L = 100 for B1 and B2 and L = 150 for B3.The network information about the geometry and the range of transmissivity values are given in Table 2  For the larger linear systems (for B1 and B2 with k = 1 and k = 2 and for B3 with k = 0, 1 and 2), we are facing with the Intel Pardiso LLT solver some problems that we are currently investigating.On the contrary, we have no problem with the Intel Pardiso LU solver.In the following Tables 3, 4 and 5, the solver information (LLT or LU) will be given for each simulation.Moreover, depending on the requirements in computational resources, the simulations have been run either on a Intel Core i7 6 cores CPU laptop (denoted by IC in the following Tables) or on a cluster node with 4 Intel Xeon E7-8890 processors and 1024 GiB of RAM (denoted by IX in the following Tables).Table 3 gives the computational time (for reading the mesh, assembling and solving the linear system) and peak memory of the NEF++ software for k = 0 for the three test cases.For the B1 test case, we provide the results obtained with the LLT and LU solvers.With LU, the RAM memory requirements are higher than with LLT, as expected.Table 3: Performance obtained with NEF++ for the HHO method with k = 0 on B1, B2 and B3.
Table 4 and Table 5 give the computational time (for reading the mesh, assembling and solving the linear system) and peak memory of the NEF++ software for k = 1 and k = 2 respectively for the three test cases.Despite B3 has fewer edges than B1, it takes more time to solve the associated linear system with a direct solver as it has more intersections (see Table 2).As shown by Tables 3, 4 and 5

Conclusion
The results in terms of computational time and accuracy we are currently obtaining with the NEF++ software are very promising to handle in a near future the millions of fractures networks provided by external industrial partners.As the RAM requirements with direct solvers are quite large, we are currently investigating the use of iterative solvers.As emphasized in this paper, the HHO method has a strong potential, also for deriving upscaled quantities owing to its high-order feature.Moreover, as HHO naturally deals with general shape elements, non-conforming discretizations at the intersections between the fractures can be naturally handled in a conforming way.Finally, as a future work, we are interested in using the dimensional templating feature offered by the DiSk++ library to solve flow in porous fractured rocks.

Figure 3 :
Figure 3: Test case B0: x-component of the equivalent permeability tensor.

Table 1 :
Table 1 presents the values of K x as the mesh is refined.Test DFN B0Equivalent permeability K x #edges MHFEM RT0 HHO, k = 0 HHO, k = 1 HHO, k = 2 case B0: 1,397 fractures, computed values of the x-component of the equivalent permeability K x with mesh refinement and different numerical methods. .

Table 2 :
Details about the three DFN test cases B1, B2 and B3.
, increasing k requires more computation times and memory as the number of dofs increases but the solutions are more accurate, as highlighted in Section 3.

Table 4 :
Performance obtained with NEF++ for the HHO method with k = 1 on B1, B2 and B3.

Table 5 :
Performance obtained with NEF++ for the HHO method with k = 2 on B1, B2 and B3.