International Conference on Adaptive Modeling and Simulation ADMOS 2009 Ph. Bouillard and P. D´ıez ©CIMNE, Barcelona, 2009

A HP FOURIER-FINITE-ELEMENT FRAMEWORK WITH MULTIPHYSICS APPLICATIONS DAVID PARDO∗ , PAWEL MATUSZYK† , MYUNG JIN NAM† ∗ Basque

Center for Applied Mathematics (BCAM), Spain e-mail: [email protected], web page: http://www.bcamath.org/pardo † Department

of Petroleum and Geosystems Engineering The University of Texas at Austin, Austin, TX, USA

Key words: Fourier Finite Element Method, Multiphysics, Goal oriented Adaptivity Abstract. We describe the design of our multiphysics high-order Fourier-Finite Element software, and we illustrate its performance by simulating petroleum engineering applications.

1

INTRODUCTION

In order to characterize a reservoir, it is customary to measure different physical phenomena both on the surface of the earth as well as at different logging positions. On the surface, oil companies typically acquire seismic and possibly marine controlled-source electromagnetic (CSEM) measurements (Fig. 1 —left— ). With these measurements, experts on the field study the expected profitability of the reservoir, possibly with the help of a numerical software based on the inversion of single-physics measurements (see, e.g. [2]). After performing this first assessment of the reservoir, various wells may be drilled into the subsurface at several locations. Subsequently, logging instruments based on different physics (such as electromagnetism, acoustics and nuclear) are introduced into each borehole, and logging measurements are acquired at various locations along the trajectory of each well. In this context, more accurate characterization of reservoirs can be obtained by employing multiphysics simulators for both direct and inverse problems. Multiphysics simulators are also needed in other industries and sciences such as the medicine (multiphysics measurements are utilized for the accurate characterization of tumors), aeronautics (multiphysics optimization is needed to design aerodynamic aircraft with low radar cross section), nano-sciences, material sciences, etc. Motivated by all these applications, in here we present a framework for solving multiphysics problems based on the use of a hp-Fourier Finite Element Method (hp-FFEM), which is the first existing multiphysics simulator based on this method. The work presented here is a continuation of [6], where we described a hp-FFEM, and [5], where we described a goal-oriented self-adaptive hp-refinement strategy. 1

David Pardo, Pawel Matuszyk, Myung Jin Nam

Figure 1: (a) Marine controlled-source electromagnetics (CSEM) scenario, and (b) logging instrument in a deviated well.

2

METHODOLOGY

Our simulator is based on a high-order Fourier-Finite Element framework suitable for the computer simulation of a large variety of 2D and 3D multiphysics problems, with special emphasis in electromagnetic (EM) applications. The resulting software incorporates a self-adaptive goal-oriented hp-refinement strategy, where element sizes h and polynomial orders of approximation p vary locally throughout the entire computational grid. The software is suitable for both shared and distributed memory parallel machines. It incorporates Lagrange—H 1 — and Nedelec—H(curl)— elements, and the implementation of Raviart-Thomas—H(div)— and L2 elements is currently under development. The framework can be employed to simulate efficiently both multiphysics as well as single-physics problems, and it works efficiently in both sequential and parallel machines. hp Fourier-Finite Element Method (hp-FFEM). This method is based on a new geometrybased formulation for simulating 3D resistivity borehole measurements employing a mix of 2D and 1D algorithms. In so doing, we utilize a 2D self-adaptive goal-oriented hp-adaptive strategy (where h indicates the element size, and p the polynomial order of approximation) combined with a Fourier series expansion in a non-orthogonal system of coordinates. This combination naturally generates a spatial domain decomposition that is used as building block for the construction of an efficient iterative solver, thereby making unnecessary the use of algebraic domain-partitioning algorithms. Moreover, the 2D self-adaptive refinement strategy enables accurate simulations of problems that include high material contrasts (occurring, for example, when simulating a metallic mandrel in an oil-based mud). Fourier Method. The Fourier transform (or Fourier series expansion in the case of finite size domains) is selected in a spatial dimension where material coefficients are as smooth as possible. Thus, the use of a high-order method is justified. We note that the Fourier dimension may not coincide with a Cartesian direction, and a possibly non-orthogonal change of coordinates may precede the selection of an appropriate dimension suitable for the Fourier method, as shown in [6]. 2

David Pardo, Pawel Matuszyk, Myung Jin Nam

A Self-Adaptive hp Goal-Oriented Algorithm. To determine an optimal distribution of element size h and polynomial order of approximation p, we employ a goal-oriented self-adaptive refinement strategy based on the iterative scheme described in [5]. At each step, given an arbitrary hp-grid, we first perform a global and uniform hp-refinement to obtain the h/2, p + 1-grid. Second, we approximate the error function in the hp-grid by evaluating the difference between the solutions associated to the hp-and h/2, p + 1-grids. If the error exceeds a user-prescribed tolerance error, then we employ the error function to guide optimal refinements over the hpgrid, and we iterate the process. Once the prescribed tolerance error has been met, we deliver the h/2, p + 1-grid as the ultimate solution of the problem. This two-dimensional refinement strategy has been proved to be efficient, robust, and highly accurate for both EM and sonic problems [3, 6]. Rham diagram. Our implementation employs multiphysics finite elements of variable order compatible with the so-called de Rham diagram [1] (Fig. 2).

R I −→ H 1   yid

∇× ∇ ∇◦ −→ H(curl) −→ H(div) −→

  yΠ

∇ R I −→ W p −→

  div yΠ

  curl yΠ

Qp

∇× −→

Vp

L2

−→

0

  yP

∇◦ −→ W p−1 −→ 0 .

Figure 2: The “de Rham” diagram is composed of two exact sequences: one at the continuous level (top), and the second one at the discrete finite element space level (bottom). Equipped with the projection based interpolation operators Π, Πcurl , and Πdiv , the “de Rham” diagram commutes.

A Perfectly Matched Layer (PML). A PML is utilized to efficiently truncate the computational domain. For details, see [4], where we demonstrate the robustness of the PML in presence of anisotropic materials at different frequencies and with high contrast in conductivity (for example, in cased wells). Solvers of Linear Equations. To ensure efficiency of the forward simulations, we need fast solvers of linear equations. We employ both direct (Gauss factorization) and iterative solvers (see Fig. 3 – left panel –). Parallel Computations. To ensure fast computations, the above method has been implemented in parallel computers for faster execution in multiple core processors as well as in parallel distributed memory machines. 3

NUMERICAL RESULTS

Left panel of Fig. 3 displays a comparison between a direct solver (MUMPS) and our iterative solver for a logging-while-drilling application. Right panel of Fig. 3 displays the solution for the marine CSEM problem described in Fig. 1 (left panel). We observe fast convergence as we increase the number of Fourier modes. 3

David Pardo, Pawel Matuszyk, Myung Jin Nam

−10

10

−11

10 Amplitude of Electric Field (V/(A m2)

Number of unknowns = 12896 * Number of Fourier Modes 250 Direct Solver Iterative Solver 200

With Oil, 0.75 Hz 1 MODE 5 MODES 9 MODES EXACT

−12

Time (sec)

10

150

−13

10

100

−14

10

50

0 0

−15

10

−16

10

5

10 15 20 Number of Fourier Modes

25

0

2000 4000 6000 8000 Horizontal Distance between TX and RX (m)

10000

Figure 3: LEFT PANEL: CPU time used by the direct and iterative solvers, respectively, as a function of the number of Fourier modes for a particular logging-while-drilling (LWD) application. For the simulations, we have employed a machine equipped with a 2.0 GHz processor and 8 GB of RAM. RIGHT PANEL: Amplitude of the electric field as a function of the horizontal distance between transmitter and receivers. Different curves indicate different numbers of Fourier modes: (a) 1 mode (dotted pink), (b) 5 modes (blue ’+’), (c) 9 modes (black circles), and (d) exact solution (red solid line). Operating frequency: 0.75 Hz.

REFERENCES [1] L. Demkowicz, P. Monk, L. Vardapetyan, and W. Rachowicz. The Rham diagram for hp finite element spaces. Computers and Mathematics with Applications, 39(7):29–38, 2000. [2] J. Gunning and M. E. Glinsky. Detection of reservoir quality using Bayesian seismic inversion. Geophysics, 72:R37–R49, 2007. [3] C. Michler, L. Demkowicz, J. Kurtz, and D. Pardo. Improving the performance of perfectly matched layers by means of hp-adaptivity. Numerical Methods for Partial Differential Equations, 23(4):832– 858, 2007. [4] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and C. Michler. PML enhanced with a self-adaptive goal-oriented hp finite-element method and applications to through-casing borehole resistivity measurements. SIAM Journal on Scientific Computing., 30:2948–2964, 2008. [5] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and L. Tabarovsky. A goal-oriented hp-adaptive finite element method with electromagnetic applications. Part I: electrostatics. International Journal for Numerical Methods in Engineering, 65:1269–1309, 2006. [6] D. Pardo, C. Torres-Verd´ın, M. J. Nam, M. Paszynski, and V. M. Calo. Fourier series expansion in a non-orthogonal system of coordinates for simulation of 3D alternating current borehole resistivity measurements. Computer Methods in Applied Mechanics and Engineering, 197:3836–3849, 2008.

4