Underground flow simulations using parallel finite element method H. Mustapha a,n , A. Ghorayeb b , K.A. Mustapha c a McGill University, Montreal, Canada, and Reservoir Engineering Research Institute, Stanford, USA b TIMC Laboratory University, Grenoble, France c Department of Mathematical Sciences, King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia article info Article history: Received 24 April 2008 Received in revised form 4 December 2008 Accepted 11 December 2008 Keywords: Finite element method Heterogeneous problem Load-balancing problem Linear flow equations Parallel equation solver abstract Advanced algorithms in parallel computing can be employed for a better understanding of groundwater flow. Water flows in a very heterogeneous geological media which contains complex structures. Decomposing these structures into approximately equivalent sub-structures for load-balancing is a major challenge. This paper proposes and analyzes a new algorithm using finite element methods to simulate parallel flow fluid in such complex media. Fully parallel software is developed and numerical examples are presented to show the efficiency and robustness of the algorithm proposed. & 2009 Elsevier Ltd. All rights reserved. 1. Introduction Natural aquifers are very heterogeneous, leading to prefer- ential flow paths and stagnant regions. Rock properties in such an aquifer change dramatically in space. Based on the spatial heterogeneity of aquifers and interrelation of effective flow properties, parallel numerical modelling is an important con- tribution to the understanding and forecasting flow behavior, and management of groundwater. This paper studies the effects of heterogeneity on the underground flow circulation through parallel simulations of flow in complex fractured media. The discrete fracture concept is used to model fractured media (Baca et al., 1984; Granet et al., 1998; Noorishad and Mehran, 1982). The fluid is assumed to only flow in the connected fractures, so- called fracture network, because the matrix is assumed to be of low permeability. This representation of fracture networks will speed up the calculation as shown (Bonnet et al., 2001; Bour and Davy, 1998; Karimi-Fard and Firoozabadi, 2003). Conventional methods used in the numerical solutions of single and multiphase flow in porous media have some limita- tions. We used the mixed finite element (MFE) method which allows for accurate calculation of fluxes; this method has been used successfully to discretize Darcy’s law (Chavent and Roberts, 1991). The mesh of a fracture network is rather complex, with interconnected 2D triangular meshes in each fracture. In this paper, we used the algorithm developed by (Mustapha, 2007; Mustapha and Mustapha, 2007) to triangulate the fractures. The discretization of the PDEs using the MFE method on the triangulation of the fractures leads to a linear system to solve with a sparse symmetric positive definite matrix. The very strong variability of hydraulic properties may lead to an ill-conditioned matrix; therefore, a direct solver is used to solve the linear system. Numerical studies for the influence of the heterogeneity need to generate a large number of realistic numerical simula- tions, leading to very large sparse linear systems. Then, two main problems have to be overcome: memory size to generate very large linear systems and CPU run time to solve a large number of linear systems. Thus, high performance computing is mandatory in this framework. In the literature, different parallel algorithms have been developed to simulate flow fluids (Wu et al., 2002a, 2002b; Zhang et al., 2003; Philip et al., 2005). These algorithms take into account the phase matrix. In this paper, the matrix is canceled out because of its low permeability. The fractures, as shown in Figs. 1 and 2, represent different sub-domain where the fluid flow through their intersections; consequently, the flow in the network can be predicted from the flow in the fractures and their intersections. For parallel simulations, the fractures can be distributed and treated separately on the processors. For numerical difficulties the fractures have not to be decomposed. The strong variation in fractures sizes constraints the load balancing on the processors. This paper presents a new algorithm to distribute heterogeneous objects (i.e. objects of different sizes) on processors, and compares it to a public algorithm (Karypis and ARTICLE IN PRESS Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/cageo Computers & Geosciences 0098-3004/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.12.018 n Corresponding author. E-mail address: hussein.mustapha@mcgill.ca (H. Mustapha). Computers & Geosciences 36 (2010) 161–166