PorePy: an open-source software for simulation of multiphysics processes in fractured porous media

Development of models and dedicated numerical methods for dynamics in fractured rocks is an active research field, with research moving towards increasingly advanced process couplings and complex fracture networks. The inclusion of coupled processes in simulation models is challenged by the high aspect ratio of the fractures, the complex geometry of fracture networks, and the crucial impact of processes that completely change characteristics on the fracture-rock interface. This paper provides a general discussion of design principles for introducing fractures in simulators, and defines a framework for integrated modeling, discretization, and computer implementation. The framework is implemented in the open-source simulation software PorePy, which can serve as a flexible prototyping tool for multiphysics problems in fractured rocks. Based on a representation of the fractures and their intersections as lower-dimensional objects, we discuss data structures for mixed-dimensional grids, formulation of multiphysics problems, and discretizations that utilize existing software. We further present a Python implementation of these concepts in the PorePy open-source software tool, which is aimed at coupled simulation of flow and transport in three-dimensional fractured reservoirs as well as deformation of fractures and the reservoir in general. We present validation by benchmarks for flow, poroelasticity, and fracture deformation in porous media. The flexibility of the framework is then illustrated by simulations of non-linearly coupled flow and transport and of injection-driven deformation of fractures. All results can be reproduced by openly available simulation scripts.publishedVersio

Une méthode mixte multi-échelles pour un simulateur de réservoir biphasé

A multiscale hybrid mixed finite element method is presented in this paper to solve two-phase flow equations on heterogeneous media under the effect of gravitational segregation. It is designed to cope with the complex geometry and inherent multiscale nature of the rocks, leading to stable and accurate multi-physics reservoir simulations. This multiscale approach makes use of coarse scale fluxes between subregions (macro domains) that allow to reduce substantially the dominant computational costs associated with the flux/pressure kernel embedded in the numerical model. As such, larger scale problems can be approximated in a reasonable computational time. Dividing the problems into macro domains leads to a hierarchy of meshes and approximation spaces, allowing the efficient use of static condensation and parallel computation strategies. The method documented in this work utilizes discretizations based on a general domain partition formed by poly-hedral subregions. The normal flux between these subregions is associated with a finite dimensional trace space. The global system to be solved for the fluxes and pressures is expressed only in terms of the trace variables and of a piecewise constant pressure associated with each subregion. The fine scale features are resolved by mixed finite element approximations using fine flux and pressure representations inside each subregion, and the trace variable (i.e. normal flux) as Neumann boundary conditions. This property implies that the flux approximation is globally H(div)-conforming, and, as in classical mixed formulations, local mass conservation is observed at the micro-scale elements inside the subregions, an essential property for flows in heterogeneous media

Discontinuous Galerkin Methods for Mass Transfer through Semi-Permeable Membranes

A discontinuous Galerkin (dG) method for the numerical solution of initial/boundary value multi-compartment partial differential equation (PDE) models, interconnected with interface conditions, is presented and analysed. The study of interface problems is motivated by models of mass transfer of solutes through semi-permeable membranes. More specifically, a model problem consisting of a system of semilinear parabolic advection-diffusion-reaction partial differential equations in each compartment, equipped with respective initial and boundary conditions, is considered. Nonlinear interface conditions modelling selective permeability, congestion and partial reflection are applied to the compartment interfaces. An interior penalty dG method is presented for this problem and it is analysed in the space-discrete setting. The a priori analysis shows that the method yields optimal a priori bounds, provided the exact solution is sufficiently smooth. Numerical experiments indicate agreement with the theoretical bounds and highlight the stability of the numerical method in the advection-dominated regime

Parallel Algorithms for the Solution of Large-Scale Fluid-Structure Interaction Problems in Hemodynamics

This thesis addresses the development and implementation of efficient and parallel algorithms for the numerical simulation of Fluid-Structure Interaction (FSI) problems in hemodynamics. Indeed, hemodynamic conditions in large arteries are significantly affected by the interaction of the pulsatile blood flow with the arterial wall. The simulation of fluid-structure interaction problems requires the approximation of a coupled system of Partial Differential Equations (PDEs) and the set up of efficient numerical solution strategies. Blood is modeled as an incompressible Newtonian fluid whose dynamics is governed by the Navier-Stokes equations. Different constituive models are used to describe the mechanical response of the arterial wall; specifically, we rely on hyperelastic isotropic and anistotropic material laws. The finite element method is used for the space discretization of both the fluid and structure problems. In particular, for the Navier-Stokes equations we consider a semi-discrete formulation based on the Variational Multiscale (VMS) method. Among a wide range of possible solution strategies for the FSI problem, here we focus on strongly coupled monolithic approaches wherein the nonlinearities are treated in a fully implicit mode. To cope with the high computational complexity of the three dimensional FSI problem, a parallel solution framework is often mandatory. To this end, we develop a new block parallel preconditioner for the coupled linearized FSI system obtained after space and time discretization. The proposed preconditioner, named FaCSI, exploits the factorized form of the FSI Jacobian matrix, the use of static condensation to formally eliminate the interface degrees of freedom of the fluid equations, and the use of a SIMPLE preconditioner for unsteady Navier-Stokes equations. In FSI problems, the different resolution requirements in the fluid and structure physical domains, as well as the presence of complex interface geometries make the use of matching fluid and structure meshes problematic. In such situations, it is much simpler to deal with discretizations that are nonconforming at the interface, provided however that the matching conditions at the interface are properly fulfilled. In this thesis we develop a novel interpolation-based method, named INTERNODES, for numerically solving partial differential equations by Galerkin methods on computational domains that are split into two (or several) subdomains featuring nonconforming interfaces. By this we mean that either a priori independent grids and/or local polynomial degrees are used to discretize each subdomain. INTERNODES can be regarded as an alternative to the mortar element method: it combines the accuracy of the latter with the easiness of implementation in a numerical code. The aforementioned techniques have been applied for the numerical simulation of large-scale fluid-structure interaction problems in the context of biomechanics. The parallel algorithms developed showed scalability up to thousands of cores utilized on high performance computing machines

Geometry–aware finite element framework for multi–physics simulations: an algorithmic and software-centric perspective

In finite element simulations, the handling of geometrical objects and their discrete representation is a critical aspect in both serial and parallel scientific software environments. The development of codes targeting such envinronments is subject to great development effort and man-hours invested. In this thesis we approach these issues from three fronts. First, stable and efficient techniques for the transfer of discrete fields between non matching volume or surface meshes are an essential ingredient for the discretization and numerical solution of coupled multi-physics and multi-scale problems. In particular L2-projections allows for the transfer of discrete fields between unstructured meshes, both in the volume and on the surface. We present an algorithm for parallelizing the assembly of the L2-transfer operator for unstructured meshes which are arbitrarily distributed among different processes. The algorithm requires no a priori information on the geometrical relationship between the different meshes. Second, the geometric representation is often a limiting factor which imposes a trade-off between how accurately the shape is described, and what methods can be employed for solving a system of differential equations. Parametric finite-elements and bijective mappings between polygons or polyhedra allow us to flexibly construct finite element discretizations with arbitrary resolutions without sacrificing the accuracy of the shape description. Such flexibility allows employing state-of-the-art techniques, such as geometric multigrid methods, on meshes with almost any shape.t, the way numerical techniques are represented in software libraries and approached from a development perspective, affect both usability and maintainability of such libraries. Completely separating the intent of high-level routines from the actual implementation and technologies allows for portable and maintainable performance. We provide an overview on current trends in the development of scientific software and showcase our open-source library utopia

Multiscale Methods for Stochastic Collocation of Mixed Finite Elements for Flow in Porous Media

This thesis contains methods for uncertainty quantification of flow in porous media through stochastic modeling. New parallel algorithms are described for both deterministic and stochastic model problems, and are shown to be computationally more efficient than existing approaches in many cases.First, we present a method that combines a mixed finite element spatial discretization with collocation in stochastic dimensions on a tensor product grid. The governing equations are based on Darcy's Law with stochastic permeability. A known covariance function is used to approximate the log permeability as a truncated Karhunen-Loeve expansion. A priori error analysis is performed and numerically verified.Second, we present a new implementation of a multiscale mortar mixed finite element method. The original algorithm uses non-overlapping domain decomposition to reformulate a fine scale problem as a coarse scale mortar interface problem. This system is then solved in parallel with an iterative method, requiring the solution to local subdomain problems on every interface iteration. Our modified implementation instead forms a Multiscale Flux Basis consisting of mortar functions that represent individual flux responses for each mortar degree of freedom, on each subdomain independently. We show this approach yields the same solution as the original method, and compare the computational workload with a balancing preconditioner.Third, we extend and combine the previous works as follows. Multiple rock types are modeled as nonstationary media with a sum of Karhunen-Loeve expansions. Very heterogeneous noise is handled via collocation on a sparse grid in high dimensions. Uncertainty quantification is parallelized by coupling a multiscale mortar mixed finite element discretization with stochastic collocation. We give three new algorithms to solve the resulting system. They use the original implementation, a deterministic Multiscale Flux Basis, and a stochastic Multiscale Flux Basis. Multiscale a priori error analysis is performed and numerically verified for single-phase flow. Fourth, we present a concurrent approach that uses the Multiscale Flux Basis as an interface preconditioner. We show the preconditioner significantly reduces the number of interface iterations, and describe how it can be used for stochastic collocation as well as two-phase flow simulations in both fully-implicit and IMPES models

Solution strategies for nonlinear conservation laws

Nonlinear conservation laws form the basis for models for a wide range of physical phenomena. Finding an optimal strategy for solving these problems can be challenging, and a good strategy for one problem may fail spectacularly for others. As different problems have different challenging features, exploiting knowledge about the problem structure is a key factor in achieving an efficient solution strategy. Most strategies found in literature for solving nonlinear problems involve a linearization step, usually using Newton's method, which replaces the original nonlinear problem by an iteration process consisting of a series of linear problems. A large effort is then spent on finding a good strategy for solving these linear problems. This involves choosing suitable preconditioners and linear solvers. This approach is in many cases a good choice and a multitude of different methods have been developed. However, the linearization step to some degree involves a loss of information about the original problem. This is not necessarily critical, but in many cases the structure of the nonlinear problem can be exploited to a larger extent than what is possible when working solely on the linearized problem. This may involve knowledge about dominating physical processes and specifically on whether a process is near equilibrium. By using nonlinear preconditioning techniques developed in recent years, certain attractive features such as automatic localization of computations to parts of the problem domain with the highest degree of nonlinearities arise. In the present work, these methods are further refined to obtain a framework for nonlinear preconditioning that also takes into account equilibrium information. This framework is developed mainly in the context of porous media, but in a general manner, allowing for application to a wide range of problems. A scalability study shows that the method is scalable for challenging two-phase flow problems. It is also demonstrated for nonlinear elasticity problems. Some models arising from nonlinear conservation laws are best solved using completely different strategies than the approach outlined above. One such example can be found in the field of surface gravity waves. For special types of nonlinear waves, such as solitary waves and undular bores, the well-known Korteweg-de Vries (KdV) equation has been shown to be a suitable model. This equation has many interesting properties not typical of nonlinear equations which may be exploited in the solver, and strategies usually reserved to linear problems may be applied. In this work includes a comparative study of two discretization methods with highly different properties for this equation

A Review of Element-Based Galerkin Methods for Numerical Weather Prediction: Finite Elements, Spectral Elements, and Discontinuous Galerkin

Numerical weather prediction (NWP) is in a period of transition. As resolutions increase, global models are moving towards fully nonhydrostatic dynamical cores, with the local and global models using the same governing equations; therefore we have reached a point where it will be necessary to use a single model for both applications. The new dynamical cores at the heart of these unified models are designed to scale efficiently on clusters with hundreds of thousands or even millions of CPU cores and GPUs. Operational and research NWP codes currently use a wide range of numerical methods: finite differences, spectral transform, finite volumes and, increasingly, finite/spectral elements and discontinuous Galerkin, which constitute element-based Galerkin (EBG) methods.Due to their important role in this transition, will EBGs be the dominant power behind NWP in the next 10 years, or will they just be one of many methods to choose from? One decade after the review of numerical methods for atmospheric modeling by Steppeler et al. (Meteorol Atmos Phys 82:287–301, 2003), this review discusses EBG methods as a viable numerical approach for the next-generation NWP models. One well-known weakness of EBG methods is the generation of unphysical oscillations in advection-dominated flows; special attention is hence devoted to dissipation-based stabilization methods. Since EBGs are geometrically flexible and allow both conforming and non-conforming meshes, as well as grid adaptivity, this review is concluded with a short overview of how mesh generation and dynamic mesh refinement are becoming as important for atmospheric modeling as they have been for engineering applications for many years.The authors would like to thank Prof. Eugenio Oñate (U. Politècnica de Catalunya) for his invitation to submit this review article. They are also thankful to Prof. Dale Durran (U. Washington), Dr. Tommaso Benacchio (Met Office), and Dr. Matias Avila (BSC-CNS) for their comments and corrections, as well as insightful discussion with Sam Watson, Consulting Software Engineer (Exa Corp.) Most of the contribution to this article by the first author stems from his Ph.D. thesis carried out at the Barcelona Supercomputing Center (BSCCNS) and Universitat Politècnica de Catalunya, Spain, supported by a BSC-CNS student grant, by Iberdrola Energías Renovables, and by grant N62909-09-1-4083 of the Office of Naval Research Global. At NPS, SM, AM, MK, and FXG were supported by the Office of Naval Research through program element PE-0602435N, the Air Force Office of Scientific Research through the Computational Mathematics program, and the National Science Foundation (Division of Mathematical Sciences) through program element 121670. The scalability studies of the atmospheric model NUMA that are presented in this paper used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. SM, MK, and AM are grateful to the National Research Council of the National Academies.Peer ReviewedPostprint (author's final draft