Research Statement
Juan Galvis
September 2011

My research interests have been related with the design and numerical analysis of effective numerical methods for flows in porous media. Applications of such a model include composite materials, hydrology, petroleum engineering, bio-fluids, and other related problems. In developing efficient numerical methods for flows, one faces a number of challenges because realistic subsurface formations vary over multiple scales and contain uncertainties. In my works, I have addressed some of these key challenges in designing efficient multiscale and domain decomposition methods for flows in porous media. My main contributions include:

• Develop and analyze robust preconditioners for flows in high-contrast multiscale media

• Develop multiscale methods for flows in high-contrast heterogeneous porous media

• Design and analyze of robust preconditioners for coupling Stokes flow and Darcy flow

• Advanced discretization schemes for coupled Stokes and Darcy flow

• Develop and analyze numerical methods for flows in random porous media

• Design and analyze robust preconditioners for Discontinuous Galerkin (DG) discretizations for flows with discontinuous coefficient.

Because of scale disparity and uncertainties, it is important to derive convergence results that are independent of the physical parameters. One of my main research goals has been to develop innovative numerical methods that converge independent of physical parameters. The research I have been involved with is strongly motivated by many fundamental issues that arise in practical applications of flows in porous media. Below I describe my research in more details.

Domain decomposition and multiscale methods for high-contrast problems

Subsurface flows are often affected by heterogeneities in a wide range of length scales. High-contrast in the media properties brings an additional small scale into the problem, expressed as the ratio between low and high conductivity values. Because of disparity of scales, it is prohibitevly expensive to do detailed fine-scale simulations. Some of the challenging issues in this are the design of reduced models on a coarse grid and robust and optimal preconditioners that converge independent of physical parameters, such as small scales and the contrast in the media properties. I have made some contributions in this direction in my papers [10,11,17,18,19]. A review of this results and some other important ideas is presented in [12].

To convey the main idea of my work, for simplicity I consider two-level domain decomposition preconditioners for flow equations

$$\left\{ \begin{array}{rcl} -\mbox{div}_{x} \left(\kappa(x)\n... ...$}\\ p(x)&=&g(x), \mbox{ on $\partial D$.} \end{array}\right.$$ (1)

Domain decomposition methods use the solutions of local problems and coarse problems in constructing preconditioners for the fine-scale system. The idea is to iterate, using local and coarse problems, until some stopping criteria is satisfied. The number of iterations required by domain decomposition preconditioners is adversely affected by the contrast in the media properties, that is, by the variation and jumps in $\kappa(x)$. Our objective is to design appropriate coarse spaces that guarantee that the preconditioners converge is independent of the contrast and thus optimal with respect to physical parameters. Previous findings in this area [27] have addressed special classes of coefficients $\kappa(x)$ (that is, special cases of distributions of heterogeneities and jumps in the coefficient).

In our papers [10,11,17,18,19] , we propose a special coarse spaces for two-level domain decomposition preconditioners that provide an optimal (w.r.t. physical parameters) preconditioners. The main idea behind the construction of these coarse spaces is the use of a local eigenvalue problem of the form $-\mbox{div}_x (k(x)\nabla_x\phi)=\widetilde{\kappa}(x)\phi$ with a modified weight function $\widetilde{\kappa}$. The weight function is computed as the pointwise energy of (multiscale) basis functions that are carefully selected. We show that the spectrum of this eigenvalue problem has a spectral gap that allows us to indentify important features of the solution and basis functions that are needed. We show that with such construction, the coarse spaces only needs to contain information and features related to long channels (high-conductivity regions that connect the boundaries of the coarse block). Coarse spaces with dimension related only to the number of channels are optimal size since it is known that one has to include the channels information in the coarse space. We prove that the convergence of two-level preconditioners is independent of physical parameters.

In [11], we discuss and analyze the coarse-scale approximation of the solution using our multiscale basis functions. In particular, we present an error analysis that shows that more basis functions are required to achieve reasonable accuracy. In [17], we extend two-level approaches to multi-level methods within the framework of spectral Algebraic Multigrid (AMG). Our main objective is to show that one can achieve optimal (in terms of physical parameters) convergence of multigrid methods if basis functions at each level are chosen appropriately. In this work, the coarse spaces are constructed hierarchically. We show that the hierarchical construction preserves important features of the solution. One of the difficulties I faced is guaranteeing that the coarse space dimension is minimal, i.e., at coarse levels, we do not represent the inclusions separately.

In [16] we show how to use our methodology to design a robust iterations for nonlinear problems (e.g., Richard's equation). In [15,14] we extended the framework introduced in [18,19,11,10,17] to other equations such as the mixed version of (1) and Brinkman. The results in [15] are general and can be applied to many equations with important physical parameters that adversely affect the performance of iterative methods. In order to apply our results to an equation with important physical parameters, we require that the problem can be formulated in a weak sense with positive bilinear forms. The dimension of the resulting coarse problem is problem dependent since the number of small eigenvalues of the resulting local generalize eigenvalue may differ for each problem. In particular, finding optimal dimension coarse space for equations such as highly anisotropic elliptic problems ([13]), elasticity, and other applications of the framework introduced in [18,19,11,10,17] is matter of current research.

Numerical approximation of random elliptic equations with lognormal distributions

Media properties contain uncertainties, especially at the fine resolution. As permeability data is also collected at finest scales such as core scales, detailed geological models are constructed that contain the information about uncertainties. At these scales, we deal with large uncertainties associated with the fine grid information and robust convergence estimates for stochastic discretization that take into account the fine-scale uncertainties are needed. In this case, it is more advantageous to work with infinite dimensional stochastic space due to a large dimension of the stochastic space.

In my work, I have considered log-normal random field so that $\xi(x,\cdot)=\log{\kappa(x,\cdot)}$ is Gaussian and its distribution is determined by its mean $E[\xi(x,\cdot)]=0$ and covariance function $C_{\xi}(x,y)=E[\xi(x,\cdot)\xi(y,\cdot)].$ This leads to the problem on a suitable probability space with infinite stochastic dimensions. A main theoretical challenge related to this log-normal coefficient is dealing with the lack of uniform ellipticity of the coefficient. One of the results in [21] is a new framework introduced to overcome this difficulty and to derive error estimates as described below.

In [35], Roman and Sarkis use the white noise framework of Hida and others to pose the problem, see [28,29,31]. The probability space is the the Schwartz space of rapidly decreasing functions $\Omega=\mathcal{S}^\prime(\mathbb{R}^d)$. They use the family $\mathcal{B}(\mathcal{S}^{\prime}(\mathbb{R}^d))$ of Borel subsets of $\mathcal{S}^{\prime}(\mathbb{R}^d)$ and the probability measure, $\mu$, is given by the Bochner-Minlos theorem. The triplet $(\mathcal{S}^{\prime}(\mathbb{R}^d),\mathcal{B}(\mathcal{S}^{\prime}(\mathbb{R}^d)),\mu)$ is the $1$-dimensional white noise probability space, and $\mu$ is called the white noise measure. The measure $\mu$ is also often called the (normalized) Gaussian measure on $\mathcal{S}^{\prime}(\mathbb{R}^d)$. Roman and Sarkis assume that the coefficient $\kappa$ is bounded away from zero. They expand the solution with respect to special basis of $L^2(\mu)$ made of multivariate (Fourier-Hermite) polynomials of random variables. This series is called the (Fourier-Hermite) chaos expansion. They proposed a finite element approximation of the problem by truncating the chaos expansion of the solutions. They provide neither a priori error estimate nor numerical experiments. In [21] we introduced adequate norms in the space $(L^2)$ and derive a priori error estimates for the Finite Element approximations proposed in [35]. See [2] and references therein. This work uses tools from the white noise calculus, also known as infinite dimensional calculus (see [28,29,31,3,34] and references therein).

In the more theoretical and very interesting submitted manuscript [26], we proved the regularity results required in [21]. In [26] we use tools from the white noise calculus and the Malliavin calculus. In particular, we prove equivalence of two important class of norms used to define regular (in $\omega$) random functions $f\in (L^2)$. The first norms are weighted chaos norms used in [21] to get a priori error estimates. The other norms are Gaussian Sobolev norms, that is, the analogous to Sobolev norms using the Gaussian measure $\mu$. This second norms is more adequate to prove regularity results, see [26].

Some numerical comparisons and possible numerical and theoretical advantages of using the framework in [21,26] are object of our current research.

Coupling fluid flow with porous media flow.

In many porous media applications, it is important to consider the coupling of free flow (described by Stokes or Navier-Stokes flow) and porous media flows. This occurs, for example, in the coupling of flow in the reservoir (described by Darcy's equations) and well bore (described by Navier-Stokes equations). The interface conditions that are often used to couple porous media and Stokes flows are Beavers-Joseph-Saffman conditions that is some type of upscaling of Stokes flow near the interface. The main challenge is to design efficient coupling of the discretization of Stokes equations and Darcy equations (e.g., multiscale discretization of Darcy described earlier). This has been a part of my research.

The purpose of my research is to analyze the coupling across an interface of fluid and porous media flows. It is considered an incompressible fluid in a region $\Omega_f$ that can flow both ways across an interface $\Gamma$ into a saturated porous medium domain $\Omega_p$. The model consists of Stokes equations in the fluid region, $\Omega_f$, and Darcy law for the filtration velocity in the porous medium region, $\Omega_p$. The transmission conditions on the interface $\Gamma$ are the Beavers-Joseph-Saffman conditions (see [1,30,36]). This model appears in several applications like well-reservoir coupling in petroleum engineering, transport of substances across groundwater and surface water, and (bio)fluid-organ interactions.

In [20,25], we consider this model. We analyze inf-sup conditions and optimal a priori error estimates associated with the continuous and discrete formulations of this Stokes-Darcy system. The continuous inf-sup analysis uses tools developed in [25] and [32]. For the discretization of this problem, Taylor-Hood and Raviart-Thomas finite elements are used for the free fluid and porous medium subdomains, respectively. Using mortar finite element analysis and appropriate scaled norms, we show that the constants that appear on the a priori error bounds do not depend on the viscosity, permeability and ratio of mesh parameters. Numerical experiments are presented to confirm the theoretical results. For standard choices of finite elements spaces and due to the small value of the permeability parameter $\kappa$ of the porous medium, the resulting discrete symmetric saddle point system is very ill conditioned. In [22,24], we design and analyze two preconditioners for the non-matching meshes case. One of the preconditioners is based on Balancing Domain Decomposition (BDD) methods and the other one based on Finite Element by Tearing and Interconnecting (FETI) methods. For both methods, we derive condition number estimates of order $C_1(1+\frac{1}{\kappa})$. In case the fluid discretization is finer than the porous side discretization, we derive a better estimate of order $C_2(\frac{\kappa+1}{\kappa+(h^p)^2} )$ for the FETI preconditioner. Here, $h^p$ is the mesh size of the porous side triangulation. The constants $C_1$ and $C_2$ are independent of the permeability $\kappa$, the fluid viscosity $\nu$, and the mesh ratio across the interface. Numerical experiments confirm the sharpness of the theoretical estimates. These two solvers are modular in the sense that they need the solution of a whole Stokes and/or Darcy problem in each iteration. This may represent a drawback of the method if these solves are not available due to time or memory budget constrains. In [23] we present some preliminary and numerical results of some (multisubdomain) FETI-DP methods for the coupling of Stokes and Darcy. These new methods allow solving the Stokes/Darcy coupled model iteratively using solution in small (either Darcy or either Stokes) subdomains. Each subdomain is a subregion of the Stokes or the Darcy type. These works use the previous multisubdomain solvers for Stokes [33] and Darcy [38]. See also references therein. The numerical analysis and other aspects of the design of methods similar to the methods we presented in [23] are under current research.

Robust solvers for discontinuous Galerkin discretization of heterogeneous elliptic equations.

The pressure equation (1) appears in the mathematical studies of the transport of pollutants in ground-water and of oil recovery processes. When $\kappa$ is piecewise constant on $D$ (or close to a piecewise constant function), a discontinuous Galerkin (DG) approximation of these elliptic problems can be used. DG methods are becoming more and more popular for the approximation of PDEs since they are well suited for dealing with regions with complex geometries or discontinuous coefficients. In my work, we have extended the design and analysis of classical domain decomposition methods to DG discretization as described below.

The region $D$ in (1) is assumed to be a geometrically conforming union of disjoint polygonal subregions $D_i$. The discontinuities of the coefficients occur across $\partial D_i$. The problem is approximated by a conforming finite element method (FEM) based on a matching triangulation inside each $D_i$, $i=1,\dots,N$. We allow these triangulations to be nonmatching across $\partial D_i$, $i=1,\dots,N$. This kind of triangulation and composite discretization is motivated, among a variety of reasons, by the regularity of the solution of the problem being discussed. The discrete problem is formulated using the symmetric DG method with interior penalty terms on $\partial D_i$, with the harmonic average values of discontinuous coefficients of the original problem on common faces of substructures. See [4,9] and references therein. In [8], we analyze Neumann-Neumann (N-N) algorithms for the resulting discrete problem. N-N methods are well known for the standard conforming and nonconforming discretizations, see [37] and references therein. This work was based on the previous work presented in a technical report by Dryja and Sarkis where they analyze several DG preconditioners of Neumann-Neumann type. We present numerical experiments that verify the theoretical estimates. In [5] and [6], we have also successfully extended these preconditioners to the Balancing Domain Decomposition (BDD) method and the Balancing Domain Decomposition with Constrains (BDDC) method, respectively. Numerical experiments are presented. In [7], we extend the BDD and BDDC algorithm for the geometrically nonconforming case. We present numerical experiments to verify the theory.

All the methods developed above assume particular distributions of the jumps of the coefficients. Our current research seeks to develop domain decomposition preconditioners for the case of more general distributions of coefficients and DG discretizations. In particular, we are studying optimal ways to apply some of the ideas introduced in [18,19,11,10,17] to this case.

Impact and Citation

According to MathSciNet and the Mathematical Reviews Database, out of my 11 entries, 7 are cited 27 times by 42 authors (by Oct 10, 2011). My MathSciNet $h$-index¹ is 3 (see here). In SciVerse Scopus, I have $h$-index 3 (see here). My Google Scholar public profile, which includes more recent posted preprint and published works, it has 134 citations with $h$-index equals to 7. The profile can be accessed here. Also, my paper [19] was listed top 20 most downloaded Articles in SIAM Multiscale Modeling and Simulation for the month of Sept, Oct, Nov, Dic (2010), Jan and Feb (2011). The same with [18] for the month Aug, Sept and Oct (2010).

Future Research

With my Phd and my experience after it, I have mastered several important modeling, numerical and mathematical tools. Also, I have become aware of many other important tools and topics in pure and applied mathematics. I am confident I can easily continue with my current research areas and current collaborations. I plan to continue developing efficient numerical methods for heterogeneous problems with multiple scales and uncertainties. As mentioned in all the topics before, there is still several research directions that I can follow up and continue my research. These directions are valid important research topics. I am also confident that I will be adding some research topics of my interest as well as starting new collaborations. Some of these topics are motivated because they are fields where I can apply the tools and areas I already know. Some others are motivated by the fact that Ill be learning and developing new tools. Possible future new research directions I consider to include into my list of research topics are: mulstiscale finite element methods for flows on rough surfaces (already started and we have results and we are preparing manuscript) ensemble-level preconditioners for stochastic problems (already started and we have nice results), Free interface problems (already started and we have some preliminary results), abstract theory for non-symmetric domain decomposition solvers, multiscale model reduction methods for heterogeneous problems, optimal preconditioners for nonlinear problems, discretization of eigenvalue problems, inverse coefficient elliptic problems, solvers and multiscale finite element methods for Stokes/Darcy coupling for the case of complex heterogeneous porous media, parabolic problems with random coefficients and temporal noise, etc.