Deep Learning for Inverting Borehole Resistivity Measurements.

Deep Learning for Inverting Borehole Resistivity Measurements Jon Ander Rivera Gonz´alez Supervised by David Pardo and Elisabete Alberdi November 2022 (cc)2022 JON ANDER RIVERA GONZALEZ (cc by 4.0) Deep Learning for Inverting Borehole Resistivity Measurements Jon Ander Rivera Gonz´alez Supervised by David Pardo and Elisabete Alberdi November 2022 This dissertation has been possible with the support of the University of the Basque Country (UPV/EHU) grant No. PIF18/017; the BCAM “Severo Ochoa” accreditation of excellence (SEV-2017-0718); the Basque Government through the BERC 2018-2021 program; and the Consolidated Research Group MATHMODE (IT1294-19;IT1456-22) given by the Department of Education. i Acknowledgements The last four years of my life have been an adventure of learning and personal growth. In this journey, I have been fortunate to have David Pardo as my supervisor. The day I met him I was surprised by the naturalness with which he treated me. From then on, I got to know the rest of his multiple virtues: work discipline (always available for any questions or queries), enthusiasm for research (always looking for new challenges and new problems to solve), and intuition–the one that impresses me the most. I will be eternally grateful to him and he knows that, wherever he is, he will always have a friend here. I would also like to thank my other supervisor Elisabete Alberdi for all the work she has done to help me. I would like to thank her for the trust she has placed in me when making me part of her projects and for all the advice she has given me during these four years. Finally, she has taught me to manage within the University and she has given me the opportunity to teach some classes and see how things look from the other side of the classroom. I am deeply grateful to Mostafa Shahriari for all his help. He is a lovely person. From the first day I came to the group, he welcomed me and had the patience to help me day after day. He taught me everything he knew and helped me build the foundation of what today is my Dissertation. I also want to thank him for the opportunity he gave me to go to Austria to collaborate with him. I learned a lot from that experience. I expect to continue seeing each other in the coming years. I would like to express my gratitude to Javier Omella for everything he has done for me. I have shared many hours of work with him and he taught me new things. He was always available to help and together we managed to get the job done. I have seen few experts advanced programmers like him and I thank him for transmitting all that knowledge to me. I wish to thank all my colleagues from the MATHMODE group and BCAM. Especially, I would like to thank Judit Mu˜noz-Matute and Ana Fernandez-Navamuel. The first one for helping me with bureaucracy questions; the second one for sharing with me her enriching experiences while developing our Dissertations. Finally, I want to thank my family and friends for their support. In particular, I want to offer my eternal gratitude to the three people that have been with me all this way long. First, to my parents Mari Carmen and Michel. Thank you for all your unconditional support and love, for teaching me the values that have made ii Acknowledgements me who I am, and for encouraging me to pursue a Ph.D degree. Last but not least, to my partner Amaia. She has consistently supported me when I needed it: she was the serenity when I was stressed and the encouragement when I was frustrated. Thank you for trusting me and for never ceasing to remind me that I was capable of doing it. Zaila bada ez da ezinezkoa, zaila bada badago lortzea. iii Abstract The Earth’s subsurface is formed by different materials, mainly porous rocks possibly containing minerals and filled with salty water and/or hydrocarbons. The formations that these materials create are often irregular, appearing geometrically abrupt forms with different properties that are mixed within the same layer. One of the main objectives in geophysics is to determine the petrophysical properties of the Earth’s subsurface. In this way, companies can discover hydrocarbon reservoirs and maximize the production, and determine optimal locations for hydrogen storage or CO2-sequestration. To achieve these goals, companies often record electromagnetic measurements using Logging While Drilling (LWD) instruments, which are able to record data while drilling. The recorded data is processed to produce a map of the Earth’s subsurface. Based on the reconstructed Earth model, the operator adjusts the well trajectory in real-time to further explore exploitation targets, including oil and gas reservoirs, and to maximize the posterior productivity of the available reserves. This real-time adjustment technique is called geosteering. Nowadays, geosteering plays an essential role in geophysics. However, it requires the capability of solving inverse problems in real time. This is challenging since inverse problems are often ill-posed. There exist multiple traditional methods to solve inverse problems, mainly, gradient-based or statistics-based methods. However, these methods have severe limitations. In particular, they often need to compute the forward problem hundreds of times for each set of measurements, which is computationally expensive in three-dimensional (3D) problems. To overcome these limitations, we propose the use of Deep Learning (DL) techniques to solve inverse problems. Although the training stage of a Deep Neural Network (DNN) may be time-consuming, after the network is properly trained, it can forecast the solution in a fraction of a second, facilitating realtime geosteering operations. In the first part of this dissertation, we investigate appropriate loss functions to train a DNN when dealing with an inverse problem. Additionally, to properly train a DNN that approximates the inverse solution, we require a large dataset containing the solution of the forward problem for many different Earth models. To create such dataset, we need to solve a Partial Differential Equation (PDE) thousands of times. Building a dataset may be time-consuming, especially for two and three-dimensional problems since solving iv Abstract PDEs using traditional methods, such as the Finite Element Method (FEM), is computationally expensive. Thus, we want to reduce the computational cost of building the database needed to train the DNN. For this, we propose the use of refined Isogeometric Analysis (rIGA) methods. In addition, we explore the possibility of using DL techniques to solve PDEs, which is the main computational bottleneck when solving inverse problems. Our main goal is to develop a fast forward simulator for solving parametric PDEs. As a first step, in this dissertation we analyze the quadrature problems that appear while solving PDEs using DNNs and propose different integration methods to overcome these limitations. v Resumen El subsuelo terrestre est´a formado por diferentes materiales, principalmente por rocas porosas que posiblemente contienen minerales y est´an rellenas de agua salada y/o hidrocarburos. Por lo general, las formaciones que crean estos materiales son irregulares y con materiales de diferentes propiedades mezclados en el mismo estrato. Uno de los principales objetivos en geof´ısica es determinar las propiedades petrof´ısicas del subsuelo de la Tierra. De este modo, las compa˜n´ıas pueden determinar la localizaci´on de las reservas de hidrocarburos para maximizar su producci´on o descubrir localizaciones ´optimas para el almacenamiento de hidr´ogeno o el dep´osito de CO2. Para este prop´osito, las compa˜n´ıas registran mediciones electromagn´eticas utilizando herramientas de Medici´on Durante Perforaci´on (MDP), las cuales son capaces de recabar datos mientras se lleva a cabo el proceso de prospecci´on. Los datos obtenidos se procesan para producir un mapa del subsuelo de la Tierra. Bas´andose en el mapa generado, el operador ajusta en tiempo real la trayectoria de la herramienta de prospecci´on para seguir explorando objetivos de explotaci´on, incluidos los yacimientos de petr´oleo y gas, y maximizar la posterior productividad de las reservas disponibles. Esta t´ecnica de ajuste en tiempo real se denomina geo-navegaci´on. Hoy en d´ıa, la geo-navegaci´on desempe˜na un papel esencial en geof´ısica. Sin embargo, requiere la resoluci´on de problemas inversos en tiempo real. Esto supone un reto, ya que los problemas inversos suelen estar mal planteados. Existen m´ultiples m´etodos tradicionales para resolver los problemas inversos, principalmente, los m´etodos basados en el gradiente o en la estad´ıstica. Sin embargo, estos m´etodos tienen graves limitaciones. En particular, a menudo necesitan calcular el problema inverso cientos de veces para cada conjunto de mediciones, lo que es computacionalmente caro en problemas tridimensionales (3D). Para superar estas limitaciones, proponemos el uso de t´ecnicas de Aprendizaje Profundo (AP) para resolver los problemas inversos. Aunque la etapa de entrenamiento de una Red Neuronal Profunda (RNP) puede requerir mucho tiempo, una vez que la red est´a correctamente entrenada puede predecir la soluci´on en una fracci´on de segundo, facilitando las operaciones de geo-navegaci´on en tiempo real. En la primera parte de esta tesis, investigamos las funciones de p´erdida apropiadas para entrenar una RNP cuando se trata de un problema inverso. vi Resumen Adem´as, para entrenar adecuadamente una RNP que se aproxime a la soluci´on inversa, necesitamos un gran conjunto de datos que contenga la soluci´on del problema directo para muchos modelos terrestres diferentes. Para crear dicho conjunto de datos, necesitamos resolver una Ecuaci´on en Derivadas Parciales (EDPs) miles de veces. La creaci´on de un conjunto de datos puede llevar mucho tiempo, especialmente para los problemas bidimensionales y tridimensionales, ya que la resoluci´on de la EDPs mediante m´etodos tradicionales, como el M´etodo de Elementos Finitos (MEF), es computacionalmente caro. Por lo tanto, queremos reducir el coste computacional de la construcci´on de la base de datos necesaria para entrenar la RNP. Para ello, proponemos el uso de m´etodos de An´alisis Isogeom´etrico refinado (AIGr). Adem´as, exploramos la posibilidad de utilizar t´ecnicas de AP para resolver EDPs, que es la limitaci´on computacional principal al resolver problemas inversos. Nuestro objetivo principal es desarrollar un simulador r´apido para resolver EDPs param´etricas. Como primer paso, en esta tesis analizamos los problemas de cuadratura que aparecen al resolver EDPs utilizando RNPs y proponemos diferentes m´etodos de integraci´on para superar estas limitaciones. vii LIST OF FIGURES 2.10 Analytical solution vs DNN predicted solution evaluated over the test dataset using the two-step based loss function. . . . . . . . . 23 2.11 Example B.1. Evolution of the different terms of the EncoderDecoder loss function given by Equation 2.25 without regularization. 26 2.12 Example B.1. Evolution of the different terms of the EncoderDecoder loss function given by Equation 2.25 with the regularization term prescribed by Equation 2.17. . . . . . . . . . . . . . . . 27 2.13 Geosignal cross-plots for the Example B.1 without regularization for the test dataset. First row: Cross-plots of type 1. Second row: Cross-plots of type 2. First column: Attenuation. Second column: Phase. ................................. 29 2.14 Cross-plots of type 4 for Example B.1 without regularization for the training dataset (first column), and with regularization for the training dataset (second column) and the test dataset (third column). First row: distance to the upper layer. Second row: distance to the lower layer. Third row: resistivity of upper layer. Fourth row: resistivity of lower layer. Fifth row: resistivity of central layer. 32 2.15 Formation of model problem I. . . . . . . . . . . . . . . . . . . . . 33 2.16 Inverted formation of model problem I using the inversion strategy of Example B.2, i.e., with input measurements corresponding to 65 logging positions per sample. . . . . . . . . . . . . . . . . . . . . . 34 2.17 Model problem I. Comparison between F˝Iand F˝Iθ˚using the inversion strategy of Example B.2, i.e., with input measurements corresponding to 65 logging positions per sample. . . . . . . . . . 35 2.18 Inverted formation of model problem I using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. . . . . . . . . . . . . . . . . . . . 37 2.19 Model problem I. Comparison between F˝Iand F˝Iθ˚without regularization using the Encoder-Decoder loss function and the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. . . . . . . . . . 38 2.20 Model problem I. Comparison between F˝Iand F˝Iθ˚using the two-step based loss function without regularization and the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. . . . . . . . . . 39 2.21 Model problem I. Comparison between F˝Iand F˝Iθ˚with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position persample. .............................. 40 xiv LIST OF FIGURES 2.22 Model problem I. Comparison between F˝Iand Fφ˚˝Iwith regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position persample. .............................. 41 2.23 Model problem I. Comparison between F˝Iand Fφ˚˝Iθ˚with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position persample. .............................. 42 2.24 Model problem 2. Comparison between actual and predicted formations with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging positionpersample. ......................... 43 2.25 Model problem 2. Comparison between F˝Iand F˝Iθ˚with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position persample. .............................. 44 2.26 Model problem 2. Comparison between F˝Iand Fφ˚˝Iwith regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position persample. .............................. 45 2.27 Model problem III, trajectory 1. Comparison between actual and predicted formations and the corresponding coaxial logs with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 47 2.28 Model problem III, Trajectory 2. Comparison between actual and predicted formations and the corresponding coaxial logs with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 48 3.1 A schematic LWD instrument with two transmitters and two receivers located symmetrically around the tool center. . . . . . . . 53 3.2 Example of the Hpcurlq ˆ H1space for a 2.5D formulation discretized by Cp´1Isogeometric Analysis (IGA) with uniform 8 ˆ8 elements in Ωx,z, polynomial degree p“4, and continuity k“3. The univariate basis functions of Hx, Hy, and Hzare shown in blue, red, and purple, respectively. Thin gray lines in the mesh skeleton denote the high-continuity element interfaces. . . . . . . . . . . . 55 xv LIST OF FIGURES 3.3 HpcurlqˆH1rIGA space in Ωx,z, associated with the 8ˆ8 domain of Figure 3.2 with p“4 and k“3, after one level of symmetric partitioning by the rIGA discretization that results in 4 ˆ4 macroelements. rIGA reduces the continuity of basis functions by k´1 degrees across the macroelement separators (the lowcontinuity bases are shown in black). Thin gray lines in the mesh skeleton denote the high-continuity element interfaces, while thick black lines illustrate the macroelement boundaries. We refer to the vertical and horizontal separators as “vs” and “hs”, respectively. . 57 3.4 A drawing of the computational domain Ωx,z and the tool trajectory. The central subgrid bounded by a magenta box is composed of a set of fine elements located in the proximity of the logging instrument. The remaining elements grow smoothly in size until reaching the boundary. . . . . . . . . . . . . . . . . . . . . . . . . 59 3.5 Numerical errors when computing attenuation ratio, A, and phase difference, P, in a homogeneous medium using IGA and rIGA discretizations, obtained by a 64 ˆ64 element mesh with different element sizes hand polynomial degrees p. ............. 61 3.6 Comparison of the decay of the numerical and analytical coaxial magnetic fields for some Fourier modes, obtained in a grid of 64ˆ64 elements with h“0.025 m and different polynomial degrees. . . . 62 3.7 Computational cost in terms of FLOPs and time for solving a 2.5D borehole resistivity problem per logging position per Fourier mode. We test rIGA discretizations with two different grids of 64ˆ64 and 128 ˆ128 elements. The computational times correspond to the use of parallel solver PARDISO using two threads. . . . . . . . . . 63 3.8 Model problem with a constant dip angle of 80˝passing through a geological fault and three different materials (well trajectory is highlighted by a red dashed line). Dimensions are in meters. . . . 65 3.9 Apparent resistivities based on the attenuation ratio, ρA, and phase difference, ρP, for the first model problem, compared with the real (exact) resistivity, ρe. We obtain the results using a rIGA discretization with 64ˆ64 elements, p“4, and 8ˆ8 macroelements. 66 3.10 Second model problem with two geological faults and inclined layers. The tool trajectory (red dashed line) has different dip angles and passes through sandstone (yellow), oil-saturated (gray), and water-saturated (green) layers. Dimensions are in meters. . . . . . 66 xvi LIST OF FIGURES 3.11 Apparent resistivities based on the attenuation ratio, ρA, and phase difference, ρP, for the second model problem, compared with the real (exact) resistivity, ρe. We employ a rIGA discretization with 64 ˆ64 elements, p“4, and 8 ˆ8 macroelements. . . . . . . . . . 67 3.12 Varying parameters at each logging position when producing the training dataset for DL inversion. . . . . . . . . . . . . . . . . . . 67 3.13 (a) Attenuation ratio vs. phase difference, and (b) apparent resistivity based on attenuation vs. apparent resistivity based on phase, obtained for the 100,000 Earth models. We use rIGA discretization with 64 ˆ64 elements, p“4, and 8 ˆ8 macroelements for generating the database. . . . . . . . . . . . . . . . . . . . . . 69 4.1 Sketch of the arquitecture of uNN ................... 77 4.2 Loss evolution of the training process for our two model problems. 78 4.3 Exact vs approximate Ritz method solutions of model problem 1 using four elements for evaluating FRpvqand a Neural Network (NN)with31weights. ........................ 79 4.4 Exact vs approximate Ritz method solutions of model problem 2 using ten elements for evaluation of FRpvqand a NN with 31 weights. 80 4.5 Exact (uexact “0) and approximated solution of a problem given by Eq. (4.27) and solved with the Least Square (LS) method. . . 80 4.6 Neural Network approximation uNN and its piecewise-linear element approximation u˚ NN,4....................... 82 4.7 Black points (dots) correspond to the original (a) training/(b) validation partitions and blue points (circles) are the points added by the refinement performed in the first and third elements. . . . . . 83 4.8 Loss evolution of the training process for our two model problems when we use a piecewise-linear approximation of the NN. . . . . . 89 4.9 Ritz method solution when we use a piecewise-linear approximation of the NN to solve the problem. . . . . . . . . . . . . . . . . 89 4.10 Loss evolution of the training process for our two model problems when using adaptive integration. . . . . . . . . . . . . . . . . . . . 90 4.11 Ritz method solution when using adaptive integration. . . . . . . 91 4.12 The solution and training information for Experiment 1 without regularization.............................. 92 4.13 The solution and training information for Experiment 1 with regularization. .............................. 93 4.14 The solution and training information for Experiment 2 without regularization.............................. 94 xvii LIST OF FIGURES 4.15 The solution and training information for Experiment 2 with regularization. .............................. 95 xviii List of Tables 2.1 Categories for geophysical variables: types Aor B. We apply a different rescaling to each of them. . . . . . . . . . . . . . . . . . 15 2.2 R2factors for cross-plots of type 1 and 2 and Examples B.1 and B.2, with and without regularization, for training and test datasets. Numbers below 0.96 are marked in boldface. . . . . . . . . . . . . 30 2.3 R2factors for cross-plots of type 3 and Examples B.1 and B.2, with and without regularization, for the test dataset. . . . . . . . 31 3.1 Computational cost for the 2.5D borehole resistivity measurements per logging position per Fourier mode. We report the solution time and FLOPs when using Cp´1IGA, rIGA with 8ˆ8 macroelements, and C0FEM with the same number of elements and polynomial degree. The computational times correspond to the use of parallel solver PARDISO using two threads. . . . . . . . . . . . . . . . . . 64 3.2 Varying parameters employed to generate the training dataset for DLinversion. ............................. 68 4.1 Loss values of the exact solution FRpuexactq, optimum piecewiselinear solution FRp˜u˚ NN,¨q(for a four and a ten equidistant element partition), and piecewise-linear solution FRpu˚ NN,¨qusing a DNN. . 88 xix 1 Introduction 1.1 Motivation and Literature Review The Earth’s subsurface is formed by different materials, mainly porous rocks containing minerals and filled with salty water and/or hydrocarbons. The formations that these materials create are often irregular, appearing abrupt forms with peaks or breaks. Furthermore, each of the several layers that compose the Earth is composed of various materials with different material properties. Figure 1.1 shows an example of a laminar subsurface formation. Figure 1.1: Earth subsurface formation with different layers. Photo taken in Sopela (Biscay, Spain). Several fields demand a map of the subsurface in order to carry out their activities, needed for: (a) minimize earthquake-induced damage, (b) enhance the production of geothermal energy, (c) store different materials such as hydrogen in subsurface reservoirs, and (d) maximize hydrocarbon recovery. In this last application, companies often record electromagnetic (EM) measurements using a Logging While Drilling (LWD) instrument. These tools incorporate different transmitters that generate an EM field. In the same way, several receivers are placed along the tool in order to receive the emitted wave after rebounding in the surroundings of the borehole. Depending on the materials and/or 1 1 Introduction the formation of the surroundings, the received waves exhibit different properties. Figure 1.2 shows an example of a conventional LWD instrument. Tx1,1 Tx1,2 Tx2,1 Tx2,2 Rx1Rx2 0.2032 m 0.8128 m,2 MHz 2.4384 m,0.25 MHz Conventional LWD Tx Rx1 Rx2 12 m,24 kHz 25 m,2 kHz Deep azimuthal Figure 1.2: Example of a conventional LWD instrument. This tool is equipped with a pair of receivers (red) and two pairs of transmitters (black). These tools have gained importance in the oil and gas industry in the last decades due to their capability to record logging data during drilling. The recorded data is processed to produce a map of the Earth’s subsurface nearby the well. Based on the reconstructed Earth model, the operator adjusts the well-trajectory in real-time to further explore exploitation targets, including oil and gas reservoirs, and to maximize the posterior productivity of the available reserves. This real-time navigation technique is called geosteering. As a consequence of the tremendous productivity increase achieved with this technique, nowadays geosteering plays an essential role in the oil and gas industry [40]. The main difficulty one faces when dealing with geosteering problems is to obtain a map of the Earth’s subsurface. We must solve the following inverse problem: given the measurements Mrecorded by the tool and the well trajectory T, we want to obtain the subsurface properties ρ. In contrast, the forward problem is the one that given the subsurface properties ρand the well trajectory T, it produces the measurements Mrecorded by the tool. Figure 1.3 presents a schematic description of the forward and inverse problems. Unfortunately, traditional inversion methods have severe limitations, which force geophysicists to continuously look for new solutions to this problem (see, e.g., [28, 41, 49, 64, 99, 128, 131, 156]). In particular, inverse problems are not well-defined, that is, there may exist multiple outputs for a given input [141, 145]. Gradient-based methods require simulating the forward problem dozens of times for each set of measurements. Moreover, these methods also estimate the derivatives of the measurements with respect to the inversion variables, which is often challenging and time consuming [141]. To alleviate the high computational costs associated with these inversion methods, simplified 1.5-dimensional (1.5D) methods are common (see, e.g., [64, 99, 133]). For the inversion of borehole resistivity measurements, an alternative is to apply statistics-based methods [53, 2 1 Introduction Subsurface properties ρ + Well trajectory T Measurements M Measurements M + Well trajectory T Subsurface properties ρ F I Forward: Inverse: Figure 1.3: Schematic description of forward and inverse problems. 85, 147]. The statistical methods also perform forward simulations hundreds of times for each set of measurements. Both gradient and statistics-based methods only evaluate the inverse operator. Thus, the entire inversion process is repeated at each new logging position. Deep Learning (DL) techniques seem apropiate to overcome the limitations of traditional methods while solving inversion problems. The large amount of research articles and industrial applications of DL algorithms in different areas – computer vision [82], speech recognition [4, 5, 154], biometrics [16], self-driving cars [54, 114], and healthcare [43, 108] to mention a few – are exponents of their high performance and capability to solve all kind of problems. In addition, in recent years there have been significant advances in the field of DL, with the appearance of Residual Neural Networks (RNNs) [59], which prevent gradient degeneration during the training stage, and Encoder-Decoder (sequenceto-sequence) Deep Neural Networks (DNNs), which have improved the DL work capability in computer vision applications [10]. Due to the high demand from industry to use DNNs, dedicated libraries and packages such as Tensorflow [82], Keras [29], and Pytorch [103] have been developed. These libraries facilitate the use of DNNs across different industrial applications [39, 66, 83, 117, 132, 144, 157]. All these advances make DNNs one of the most powerful and fast-growing Artificial Intelligence (AI) tools presently. The first main contribution of this dissertation is to design a fast inversion method using DL techniques to solve borehole measurement problems that allows the application of geosteering techniques. However, DNNs also face important challenges when applied to the inversion of borehole resistivity problems. In particular, the training stage can be time3 1 Introduction consuming. However, this is an offline cost incurred during the training stage. Then, after the network is properly trained, it can forecast a solution in a fraction of a second [131]. This feature allows real-time inversion, which facilitates geosteering operations. Another limitation of DNN is that they require a large dataset (also known as ground truth). In our case, it consists of the solution of the forward problem for different Earth models [60, 130, 131]. To generate the database for DL inversion, we must solve the forward problem. Our forward problem – simulation of borehole measurements – is governed by Partial Differential Equations (PDEs). In our case, we consider resistivity measurements governed by a set of four time-dependent first-order PDE named Maxwell’s equations [45]. Here, knowing some electrical properties (i.e., electrical conductivities of the subsurface materials), we can obtain the corresponding electric and magnetic fields (i.e., recorded measurements). We solve the forward problem using numerical simulation methods such as the Finite Element Method (FEM) [6, 11, 58, 69, 98, 115, 133] or the Finite Difference Method (FDM) [35, 36, 77, 140]. Moreover, we need to optimally sample the parameter space describing relevant Earth models. This process may be timeconsuming, especially for two and three-dimensional problems. In those cases, it is common to reduce the Earth model dimensionality to two or one spatial dimensions using a Fourier or a Hankel transform. These transformations lead to the so-called 2.5D [3, 48, 92, 101, 135] and 1.5D [12, 99, 129] formulations, respectively. In particular, 1.5D simulations are inaccurate when dealing with geological faults. Galerkin methods are effective for simulating well-logging problems (see, e.g., [24, 27, 84, 95, 97, 120, 146]). Isogeometric Analysis (IGA), introduced by [61], is a widely used Galerkin method for solving PDEs. IGA has been successfully employed in various EM [20, 21, 93, 137, 138] and geotechnical [56, 134] applications. IGA uses spline basis functions introduced in Computer-Aided Design (CAD) as basis functions of FEM. These basis functions exhibit high continuity (up to Cp´1, being pthe polynomial order of spline bases) across the element interfaces. When comparing IGA and FEM, the former provides smoother solutions for wave propagation problems with a lower number of unknowns [31, 61]. However, in contrast to the minimal interconnection of elements in FEM, high-continuity IGA discretizations strengthen the interconnection between elements, leading to an increase of the cost of matrix LU factorization per degree of freedom when using sparse direct solvers [30]. In order to avoid this degradation and also benefit from the recursive partitioning capability of multifrontal direct solvers, [47] developed a new method called refined Isogeometric Analysis (rIGA). This discretization technique conserves desirable properties of high-continuity IGA discretizations, while it partitions the computational domain into blocks of macroelements 4 2 Solving Inverse Problems using Deep Learning receiver is significantly larger than that of the previously considered LWD instrument. It also employs tilted receivers that are sensitive to the presence of bed boundaries. We record several measurements with this logging instrument: (a) the attenuation and phase differences, denoted as deep coaxial, computed using Equation (2.3) with H2 zz “1, and (b) the attenuation and phase differences of a directional measurement expressed as: Geosignal “ln Hzz ´Hzx Hzz `Hzx “ln |Hzz ´Hzx | |Hzz `Hzx | loooooooomoooooooon ˆ20 logpeq“:attenuation pdBq `ipphpHzz ´Hzxq´phpHzz `Hzxqq looooooooooooooooooooomooooooooooooooooooooon ˆ180 π“:phase difference (degree) . (2.4) These measurements exhibit a discontinuity as a function of the dip angle at 90 degrees. Indeed, such discontinuity is essential in the measurements if one wants to discern between top and bottom of the logging instrument (see Figure 2.4). 𝑇 𝑥 𝑅𝑥𝑇 𝑥 𝑅𝑥 𝑇 𝑥 𝑅𝑥𝑇 𝑥 𝑅𝑥 100 𝛺⋅𝑚 1𝛺⋅𝑚 100 𝛺⋅𝑚 𝐷 𝐶 𝐵 𝐴 Figure 2.4: Illustration with four logging trajectories. By symmetry, measurements recorded with trajectories A and D are identical. The same occurs with trajectories B and C. If these measurements are continuous with respect to the dip angle, then they coincide at 90 degrees, which disables the possibility of identifying if a nearby bed boundary is located above or below the logging instrument. For our borehole resistivity applications, we consider a zero-thickness borehole embedded in a three-layer medium (see Figure 2.5). A common practice in the field is to characterize this medium with seven parameters, as described in Figure 2.5. In this work, to simplify the problem, we consider only five of them by restricting the search to isotropic formations (ρv“ρh) with zero dip angle (β“ 0), as illustrated in Figure 2.6. Thus, np“5. 11 2 Solving Inverse Problems using Deep Learning Borehole 𝑑𝑢 𝑑𝑙 𝛽 𝜌𝑙 𝜌ℎ 𝜌𝑣 𝜌𝑢 Figure 2.5: Well trajectory in a 1D medium. The black circle indicates the last trajectory position. ρhand ρvare the horizontal and vertical resistivities of the host layer corresponding to the final logging position, respectively. ρuand ρlare the resistivity values of the upper and lower layers to the host layer, respectively. duand dlshow the distance from the final logging position to the upper and lower bed boundaries, respectively. t 𝑑𝑢 𝑑𝑙 𝜌𝑢 𝜌𝑙 𝜌ℎ (a) Example B.1: trajectory with 1 logging positions 𝐭 𝑑𝑢 𝑑𝑙 𝜌𝑢 𝜌𝑙 𝜌ℎ (b) Example B.2: trajectory with 65 logging positions Figure 2.6: Model problems corresponding to examples B.1 and B.2, respectively. In this example, we consider two cases (see Figure 2.6) according with the different numbers of logging positions we consider per data sample. 2.1.5.1 Example B.1: One Logging Position In this case, each trajectory consists of a single logging position. Therefore, for each sample, we record six real numbers (three attenuations and three phases), i.e., nm“6. At each logging position, the trajectory is described by one number: the trajectory dip angle. Thus, nt“1. 12 2 Solving Inverse Problems using Deep Learning 2.1.5.2 Example B.2: Sixty-Five Logging Positions In this case, the logging trajectory of each sample is formed by 65 logging positions with a logging step size of 0.3048 m(see [130, 131] for further details). Thus, for each Earth model p, we parametrize mwith 6 ˆ65 “390 real numbers (nm“390). For this example, we assume that the variation of the dip angle at a given logging position with respect to the previous one is constant. We denote that constant dip angle variation as αv. Then, at the i-th logging position, the trajectory dip angle is αi“αini ` pi´1qαv, where αini is the initial dip angle. Hence, we have nt“2. 2.2 Data Space and Ground Truth In this work, we employ a DNN to approximate the discrete inverse operator I. Given a supervised database of n-pairs pmi,Ipti,miqq,i“1, ..., n, the DNN builds an approximation of the unknown function I. This section describes the construction of the supervised database. We first select the number of samples, n, and two subspaces of Rnpand Rnt, respectively. Then, we select the nsamples in those subspaces, namely, ppt1,p1q, ..., ptn,pnqq. To each of these samples, we apply the operator F. That is, we compute pFpt1,p1q, ..., Fptn,pnqq. Finally, the n-pairs pmi,Ipti,miqq :“ pFpti,piq,piq,i“1, ..., n form our supervised database. We denote by TPRntˆnto the set of all trajectory samples pt1, ..., tnq. In other words, Tis a matrix with tibeing its i-th column. Similarly, we define M“ pm1, ..., mnq P Rnmˆnand P“ pp1, ..., pnq P Rnpˆn. Example A: Simple model problem with known analytical solution. We select n“103uniformly spaced samples within the subspace r´33,33s Ă R. Example B: Inversion of borehole resistivity measurements. We select n“ 106. Then, for the five parameters described in Section 2.1.5, we select random samples of the following rescaled variables over the corresponding intervals forming a subspace of R5: logpρlq,logpρuq,logpρhq P r0,3s logpdlq,logpduq P r´2,1s.(2.5) We consider arbitrary high-angle trajectories. For each model problem, we randomly select the trajectory parameters within the following intervals: αini P r830,970s αvP r´0.0450,0.0450s (only for Example B.2).(2.6) 13 2 Solving Inverse Problems using Deep Learning 2.3 Data Preprocessing Notation. For each output parameter of Fand I, we denote by x“ px1, ..., xnq the n-samples associated with that parameter. These xiare real scalar values for i“1, ..., n. For example, in the borehole resistivity example, each variable x contains nsamples of each particular geophysical quantity such as resistivities, distances, or given measurements (attenuations, phases, etc.). Each dimension corresponds to a particular value (sample) of that variable, for example, the geosignal attenuation recorded at a specific logging position. From the algebraic point of view, the variable xdenotes a row of either matrix Mor P. Data preprocessing algorithm. This algorithm consists of three steps. 1. Logarithmic change of coordinates. We introduce the following change of variables: Rlnpxq:“ pln x1, ..., ln xnq.(2.7) For some geophysical variables (e.g., resistivity), this change of variables ensures that equal-size relative errors correspond to similar-size absolute errors. Thus, this change of variables allows us to perform local (within a variable) comparisons. 2. Remove outlier samples. In practice, often outlier measurements are present in the sample database. These outliers appear due to measurement error or the physics of the problem. For example, in borehole resistivity measurements, some apparent resistivity measurements approach infinity, producing “horns” in the logs. When outlier measurements exists in any particular variable of the i-th sample xi, then the entire sample should be removed. Otherwise, outlier measurements affect the entire minimization problem, leading to poor numerical results. The removal process may be automated using statistical indicators, or decided by the user based on a priori physical knowledge about the problem. We follow this second approach in this work. 3. Linear change of coordinates. We now introduce a linear rescaling mapping into the interval r0.5,1.5s. We select this interval since it has unit length and the mean of a normal (or a uniform) distribution variable xis equal to one. Let xmin :“minixi,xmax :“maxixi. We define Rlinpxq:“ˆx1´xmin xmax ´xmin `0.5, ..., xn´xmin xmax ´xmin `0.5˙,(2.8) where the limits xmin and xmax are fixed for all possible approximations xapp. This change of variables allows us to perform a global comparison between 14 2 Solving Inverse Problems using Deep Learning errors corresponding to different variables since they all take values over the same interval. Remark: xmin and xmax could also be selected based on the physically valid interval of each particular variable rather than on the training samples. Variables classification. We categorize each input and output geophysical variable xinto two types: either linear (A) or log-linear (B). When necessary, we shall indicate that a particular variable belongs to a specific category by adding the corresponding symbol as subindex of the variable, e.g., xA. Table 2.1 describes the domain of those variables as well as the rescaling employed for each of them. Variables of type Aonly require a global rescaling while those of type Brequire both a local and a global change of variables. Geophysical Variables Category Domain Rescaling Angles, attenuations, ARnRlinpxq phases, and geosignals Apparent resistivities, Bpa, 8qnRlinpRlnpxqq resistivities, and distances aą0 Table 2.1: Categories for geophysical variables: types Aor B. We apply a different rescaling to each of them. For simplicity, we denote by Rthe result of the above rescalings, i.e., RpxAq:“ RlinpxAq, and RpxBq:“RlinpRlnpxBqq. In general, given a variable x(of category Aor B), we represent xR:“Rpxq. Given a matrix XPRnxˆn, we abuse notation and denote by XR:“RpXq P Rnxˆnto the matrix that results from applying operator Rrow-wise. Remark: Substituting in Equation 2.7 the natural logarithm by the base ten logarithm does not affect the definition of R. Results are identical. 2.4 Norms and Errors We first introduce both the vector and the matrix norms that we use during the training process. Norms. We introduce a norm ||¨||Xassociated with the variable x. In general, we employ the l1or l2vector norms and, for matrices, the l1and Frobenius norms. 15 2 Solving Inverse Problems using Deep Learning Absolute and relative errors. Let xapp “ pxapp 1, ..., xapp nqbe an approximation of x. We define the absolute error Aebetween xapp and xin the ||¨||Xnorm as AX epxapp,xq:“ ||xapp ´x||X.(2.9) This error measure has limited use since it is challenging to select an absolute error threshold that distinguishes between a good and a bad quality approximation. To overcome this issue, practitioners often employ relative errors. We define the relative error Rein percent between xapp and xin the ||¨||Xnorm as: RX epxapp,xq:“100||xapp ´x||X ||x||X .(2.10) Error control. For a variable xand its approximation xapp, we want to control the relative error of the rescaled variable, that is: RX epxapp R,xRq.(2.11) The value B“ ||xR||Xis expected to be similar for all variables x. Thus: ÿ xR AX epxapp R,xRq “ ÿ xR||xapp R´xR||X«Bÿ xR ||xapp R´xR||X ||xR||X“B 100 ÿ xR RX epxapp R,xRq. (2.12) Therefore, the minimum of the first and last terms of the above equation coincide. 2.5 DNN Architectures To approximate the forward and inverse problems, we use DNN architectures based on residual-type blocks [59, 110] with convolutional operators [60, 67, 74, 151]. This work does not discuss optimal data sampling techniques nor the decision-making for the optimal selection of DNN architectures [89, 109]. In the following, we first define the main operators of our DNN architectures, followed by a description of the forward and inverse DNN architectures. We denote by Nto our nonlinear activation function. In our case, we employ the rectified linear unit (ReLU), defined component-wise for each entry x as maxp0, xq[60]. We now introduce a 1D convolutional operator Cc,k ψ, where cis the filter size (output dimensionality), kthe kernel size, and ψthe weights [59, 60]. Then, we define the following block: Bc,k ψ:“´N˝Cc,k ψ1˝N˝Cc,k ψ2`Cc,k ψ3¯,(2.13) where now ψ“ pψ1, ψ2, ψ3qare all weights associated to block Bc,k ψ. 16 2 Solving Inverse Problems using Deep Learning 2.5.1 Forward Problem DNN Architecture Each input sample has dimension np`nt, and contains the variables representing the material properties and the trajectory. We define our DNN architecture as: FR,φ :“N˝Cc6,1 φ6˝L˝U˝Bc5,3 φ5˝U˝Bc4,3 φ4˝¨¨¨˝U˝Bc1,3 φ1,(2.14) where •Uis a 1D upsampling operator with upsampling factor equal to two (using the TensorFlow routine upsampling1D [2, 102, 142]). •Lis a bilinear resampling operator with resampling factor equal to the number of logging positions [102, 142], i.e., 1 for Example B.1 and 65 for Example B.2. •ci:“40i, for i“1,¨¨¨ ,5 and c6“n1 m“6, where n1 mis the number of evaluated measurements per logging position. •φ“ tφi:i“1,¨¨¨ ,6uis a set of all weights associated to the forward DNN architecture. Lexpands (in case of 65 logging positions) or shrinks (in case of one logging position) its input dimension. The output of the mentioned bilinear resampling is a matrix in which its first dimension is equal to the number of logging positions [102, 142]. All the resampling operators considered in the Equation 2.14 raise/shrink the dimension of their input gradually to avoid missing information due to a sudden dimension change. The output is a matrix of dimension pnl, n1 mq, where nlis the number of logging positions. 2.5.2 Inverse Problem DNN Architecture The input of the DNN is a matrix of dimension pnl, n1 m`2q, where nlis the number of logging positions. The first two columns of the aforementioned matrix are the sine and cosine of the trajectory dip angle at each logging position. Analogously to the forward problem, we consider the following architecture: IR,θ :“N˝Dnp θ7˝S˝Bc6,3 θ6˝Bc5,3 θ5˝Bc4,3 θ4˝¨¨¨˝Bc1,3 θ1,(2.15) •Dn θis a fully-connected layer with nbeing its number of units and θits weights [2, 60]. •Sis a flattening layer that receives a 2D matrix and outputs a 1D vector [2]. 17 2 Solving Inverse Problems using Deep Learning •ci:“40i, for i“1,¨¨¨ ,5. •θ“ tθi:i“1,¨¨¨ ,7uis a set of all the weights associated to each block and layer Dnp θ7performs the ultimate feature extraction and down-sampling. The output of this DNN is a vector consisting of material properties. 2.6 Loss Function In this section, we consider a set of weights θPΘ and a function IR,θ that depends upon the selected DNN architecture. Then, we introduce a loss function LpIR,θq. We define the minimizer of the loss function over all possible weight sets θas: IR,θ˚:“arg min θPΘLpIR,θq.(2.16) Function Iθ˚:“R´1˝IR,θ˚˝Ris the final DNN approximation of I. In the following, we analyze the advantages and limitations associated with the use of different loss functions. 2.6.1 Data Misfit A simple loss function based on the data misfit is given by: LpIR,θq:“ ||IR,θpTR,MRq´PR||P.(2.17) In the above equation, symbol ||¨||Pindicates l1or Frobenius norms introduced in Section 2.4. Example A: Model problem with known analytical solution. In this example, np“1. Thus, matrix norms reduce to vector norms. Figure 2.7 illustrates the results we obtain using the l1and l2norms, respectively. These disappointing results are expected. In the case of l2-norm, for a sufficiently flexible DNN architecture the exact solution is Iθ˚“0. We want to minimize ÿ iPIpIR,θpmiq´piq2,(2.18) where I“ t1, ..., nudenotes the training dataset. For every sample of the form pmi,?miq, there exists another one pmi,´?miq, which is satisfied in our dataset 18 2 Solving Inverse Problems using Deep Learning by construction (see Section 2.2). Then, for a specific sample mi, the solution that minimizes the loss must satisfy pIR,θpmiq´?miq2`pIR,θpmiq´p´?miqq2.(2.19) Taking the derivative of Eq. (2.19) with respect to IR,θpmiqand equaling it to zero, we obtain 4¨IR,θpmiq “ 0.(2.20) Thus, for any sample mi, the function is minimized when the approximated value is IR,θpmiq “ 0. We thus conclude that for l2-norm, the approximated solution must be Iθ˚“0. In the case of l1-norm, any solution in between the two square root branches is valid. We want to minimize ÿ iPI|IR,θpmiq´pi|,(2.21) where I“ t1, ..., nudenotes the training dataset. Then, for a specific sample mi the solution that minimizes the loss must satisfy |IR,θpmiq´?mi|`|IR,θpmiq´p´?miq|.(2.22) By analyzing each possible case, we can express Eq. (2.22) as follows $ & % ´2¨IR,θpmiq,if IR,θpmiqă´?mi, 2?mi,if ´?miďIR,θpmiq ď ?mi, 2¨IR,θpmiq,if IR,θpmiq ą ?mi. (2.23) We see that the loss function for the exact solution attains its minimum at every point Iθ˚pmiq P r´?mi,?mis. Moreover, our numerical solutions in Figure 2.7 confirm these simple mathematical observations. Thus, the data misfit loss function is unsuitable for inversion purposes. 2.6.2 Misfit of the Measurements To overcome the aforementioned limitation, we consider the following loss function that measures the misfit of the measurements (see [68]): LpIR,θq:“ }pFR˝IR,θqpTR,MRq´MR}M,(2.24) where FR:“R˝F˝R´1, and || ¨ ||Mindicates a matrix norm of the type introduced in Section 2.4. 19 2 Solving Inverse Problems using Deep Learning 0 500 1,000 −33 0 33 Predicted Real 𝑚 𝜃∗(𝑚) (a) }¨}1-norm 0 500 1,000 −33 0 33 Predicted Real 𝑚 𝜃∗(𝑚) (b) }¨}2-norm Figure 2.7: Analytical solution vs DNN predicted solution evaluated over the test dataset using the loss function based on the data misfit. 0 500 1,000 −33 0 33 Predicted Real 𝑚 𝜃∗(𝑚) (a) }¨}1-norm 0 500 1,000 −33 0 33 Predicted Real 𝑚 𝜃∗(𝑚) (b) }¨}2-norm Figure 2.8: Analytical solution vs DNN predicted solution evaluated over the test dataset using the loss function based on the measurements misfit. 20 2 Solving Inverse Problems using Deep Learning 0 400 800 10−2 Epoch nr. Loss training validation (a) }FR,φpTR,PRq ´ MR}M 0 400 800 10−2 Epoch nr. Loss training validation (b) }pFR,φ ˝IR,θqpTR,MRq ´ MR}M 0 400 800 10−1 Epoch nr. Loss training validation (c) ||IR,θpTR,MRq ´ PR||P 0 400 800 10−1 Epoch nr. Loss training validation (d) Total Loss Figure 2.12: Example B.1. Evolution of the different terms of the EncoderDecoder loss function given by Equation 2.25 with the regularization term prescribed by Equation 2.17. 27 2 Solving Inverse Problems using Deep Learning we can safely omit these cross-plots. Otherwise, cross-plots display interesting information beyond what R2provides. The proper interpretation of the cross-plots (or alternatively, R2factors) is of utmost importance. Cross-plots of type 1 (Equation 2.281) indicate how well the forward function is approximated over the given dataset. The cross-plots of type 2 (Equation 2.282) display how well the composition of the predicted forward and inverse mappings approximate the identity. These two types of cross-plots often deliver high R2factors, since the corresponding approximations are directly built into the Encoder-Decoder loss function given by Equation 2.25. Table 2.2 confirms those theoretical predictions for the most part. An in-depth inspection of Table 2.2 reveals that for the the geosignal measurements (both attenuation and phase) corresponding to the Example B.1 without regularization, the cross-plots 2 exhibit significantly better R2factors than those corresponding to the cross-plots 1. Figure 2.13 shows the corresponding crossplots. The anti-diagonal grey line shown in cross-plots of type 1 corresponds to dip angles of the logging instrument that are close to 90 degrees. At that angle, the geosignal is discontinuous. Thus, it is not properly approximated via DL algorithms, which approximate continuous functions. Cross-plots of type 2 seem to fix that issue by delivering higher R2factors and apparently nicer figures. However, they amplify the problem. In reality, the DL approximation of the inverse operator is inverting an incorrect forward approximation. Numerical results below illustrate this problem. Obtaining high R2factors associated to cross-plots of type 3 (Equation 2.283) is a challenging task as we discuss in Remark A of Section 2.6. Equation 2.27 shows a simple example in which cross-plots of type 1 and 2 deliver perfect R2marks and results, while cross-plots of type 3 are disastrous. This is also the situation that occurs in Example B.2. (see Table 2.3). While the original training dataset is based on 1D Earth models, the one obtained after the predicted DNN inversion is a piecewise 1D Earth model, for which Fφ˚is untrained for. When this occurs, the training database should be upgraded, either by increasing the space of the data samples or by selecting a different parameterization (e.g., measurements) for each sample. In our case, we choose to parametrize each sample independently (the later stategy) and we move to Example B.1. Table 2.3 shows mixed results for the Example B.1. Results without regularization are unremarkable with the geosignal forecasts showing poor results. The DNN inverse approximation accurately inverts for the outcome predicted by the DNN forward approximation. Nevertheless, since the DNN predicts solutions far from the true forward function, the predictions are poor. Again, this poor forecasting occurs because the DNN inverse approximation encounters subsurface models for which the forward DNN approximation is untrained. As a result, 28 2 Solving Inverse Problems using Deep Learning 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Ground Truth 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Predicted r2= 0.9531 Atten-Geosignal −1 0 1 −1 0 1 Ground truth Predicted value 0.75 0.50 0.25 0.00 0.25 0.50 0.75 Ground Truth 0.75 0.50 0.25 0.00 0.25 0.50 0.75 Predicted r2= 0.9487 Phase-Geosignal −1−0.5 0 0.5 1 −1 −0.5 0 0.5 1 Ground truth Predicted value 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Ground Truth 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Predicted r2= 0.9999 Atten-Geosignal_FI −1 0 1 −1 0 1 Ground truth Predicted value 0.75 0.50 0.25 0.00 0.25 0.50 0.75 Ground Truth 0.75 0.50 0.25 0.00 0.25 0.50 0.75 Predicted r2= 0.9999 Phase-Geosignal_FI −1−0.5 0 0.5 1 −1 −0.5 0 0.5 1 Ground truth Predicted value Figure 2.13: Geosignal cross-plots for the Example B.1 without regularization for the test dataset. First row: Cross-plots of type 1. Second row: Cross-plots of type 2. First column: Attenuation. Second column: Phase. 29 2 Solving Inverse Problems using Deep Learning Cross-plots 1 Atten. Atten. Atten. Phase Phase Phase R2factors LWD Deep Deep LWD Deep Deep Coaxial Coaxial Geosignal Coaxial Coaxial Geosignal Example B.1 Training 0.9997 0.9992 0.9509 0.9996 0.9994 0.9468 Test 0.9995 0.9984 0.9531 0.9990 0.9991 0.9487 Without Reg. Example B.1 Training 0.9998 0.9998 0.9897 0.9998 0.9998 0.9893 Test 0.9998 0.9998 0.9893 0.9998 0.9998 0.9890 With Reg. Example B.2 Training 0.9959 0.9975 0.9872 0.9954 0.9980 0.9853 Test 0.9924 0.9960 0.9775 0.9920 0.9974 0.9765 Without Reg. Cross-plots 2 Atten. Atten. Atten. Phase Phase Phase R2factors LWD Deep Deep LWD Deep Deep Coaxial Coaxial Geosignal Coaxial Coaxial Geosignal Example B.1 Training 0.9997 0.9995 0.9998 0.9999 0.9996 0.9999 Test 0.9997 0.9994 0.9999 0.9999 0.9996 0.9999 Without Reg. Example B.1 Training 0.9971 0.9980 0.9779 0.9970 0.9979 0.9798 Test 0.9970 0.9979 0.9785 0.9970 0.9978 0.9803 With Reg. Example B.2 Training 0.9931 0.9958 0.9800 0.9933 0.9967 0.9821 Test 0.9890 0.9930 0.9701 0.9881 0.9944 0.9720 Without Reg. Table 2.2: R2factors for cross-plots of type 1 and 2 and Examples B.1 and B.2, with and without regularization, for training and test datasets. Numbers below 0.96 are marked in boldface. 30 2 Solving Inverse Problems using Deep Learning both the forward and inverse DDN approximations depart strongly from the true solutions. In other words, the inverse can only comply with their composition to be close to the identity, which is not robust to deliver accurate and physically relevant approximations. Cross-plots 3 Atten. Atten. Atten. Phase Phase Phase R2factors LWD Deep Deep LWD Deep Deep Coaxial Coaxial Geosignal Coaxial Coaxial Geosignal Example B.1 Without Reg. 0.9468 0.7406 0.0013 0.9383 0.9116 0.0167 With Reg. 0.9971 0.9979 0.9807 0.9969 0.9979 0.9856 Example B.2 Without Reg. 0.5721 0.8383 0.0253 0.4546 0.8611 0.0284 With Reg. 0.9010 0.9701 0.5901 0.8621 0.9618 0.5877 Table 2.3: R2factors for cross-plots of type 3 and Examples B.1 and B.2, with and without regularization, for the test dataset. To partially alleviate the above problem, we envision three possible solutions. First, we can increase the training dataset. This option is time-consuming and often impossible to achieve in practice. For example, herein, we already employ 1,000,000 samples. Second, we can include regularization. Results with regularization are of high quality (see Table 2.3). However, the regularization term may hide alternative physical solutions of the inverse problem. Thus, the regularization diminishes the ability to perform uncertainty quantification. Similarly, it may induce on the user excessive confidence in the results. A third option is to consider the two-step based loss function given by Equation 2.26. Following this approach, we first adjust the forward DNN approximation before training the DNN inverse approximation. Fixing the forward DNN often provides a proper forecast even in areas with a lower rate of training samples before producing a DNN approximation that approximates the inverse of the DNN forward approximation. Following this two-step approach without regularization, we obtain high R2factors for cross-plots of type 3: above 0.95 for the geosignal attenuation and phase, and above 0.99 for the remaining measurements. Finally, the R2factors for the cross-plots of type 4 do not reflect on the accuracy of the DNN algorithm, but rather on the nature of the inverse problem at hand. Low R2factors indicate there exist multiple solutions. A regularization term (e.g., Equation 2.17) increases the R2indicator. Figure 2.14 clearly illustrates this fact. However, it is misleading to conclude that results without regularization 31 2 Solving Inverse Problems using Deep Learning are always worse. They may simply exhibit a different (but still valid) solution of the inverse problem. 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.0227 d_u −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.7488 d_u −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.7488 d_u −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.0000 d_l −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.7466 d_l −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Ground Truth 2.0 1.5 1.0 0.5 0.0 0.5 1.0 Predicted r2= 0.7466 d_l −2−1 0 1 −2 −1 0 1 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.1876 rho_u 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.9732 rho_u 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.9706 rho_u 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.0643 rho_l 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.9742 rho_l 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.9717 rho_l 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.1508 rho_h 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.8950 rho_h 0123 0 1 2 3 Ground truth Predicted value 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Ground Truth 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Predicted r2= 0.8910 rho_h 0123 0 1 2 3 Ground truth Predicted value Figure 2.14: Cross-plots of type 4 for Example B.1 without regularization for the training dataset (first column), and with regularization for the training dataset (second column) and the test dataset (third column). First row: distance to the upper layer. Second row: distance to the lower layer. Third row: resistivity of upper layer. Fourth row: resistivity of lower layer. Fifth row: resistivity of central layer. 32 2 Solving Inverse Problems using Deep Learning 2.8.3 Inversion of Realistic Synthetic Models We now consider three realistic synthetic examples to assess the performance of the inversion process. In terms of log accuracy, we observe qualitatively similar results for the attenuation and phase logs. Thus, in the following we only display the attenuation logs and omit the phase logs. 2.8.3.1 Model Problem I Figure 2.15 describes a well trajectory in a synthetic model problem. The model has a resistive layer with a water-bearing layer underneath, and exhibits two geological faults. 0 50 100 150 200 250 300 350 400 450 500 45 50 55 60 HD (𝑚) TVD (𝑚) Figure 2.15: Formation of model problem I. For the DNNs produced with the Example B.2 (with input measurements corresponding to 65 logging positions per sample), Figure 2.16 shows the corresponding inverted models using the Encoder-Decoder DNN with and without regularization. Results show inaccurate inversion results, specially for the case without regularization. Moreover, the predicted logs are far from the true logs, as Figure 2.17, and as expected from cross-plots 3 (see Table 2.3). The DNN inversion results are piecewise 1D models. However, the DNN approximation only trains with 1D models, not for piecewise 1D models, which explains the poor approximations they deliver (see Remark A on Section 2.6). In the remainder of this section, we restrict to DNNs produced with Example B.1. That is, we parametrize all observations at one location using information from that location alone. Figure 2.18 shows the corresponding inverted models. For the case of the Encoder-Decoder loss function without regularization, we observe in Figure 2.18a an inverted model that is completely different from the original one. The corresponding logs (see Figure 2.19) are also inaccurate, as anticipated by the cross-plots results of type 3 shown in the previous subsection. When considering the two-step based loss function without regularization, the 33 2 Solving Inverse Problems using Deep Learning 0 50 100 150 200 250 300 350 400 450 500 45 50 55 60 HD (𝑚) TVD (𝑚) (a) Without regularization 0 50 100 150 200 250 300 350 400 450 500 45 50 55 60 HD (𝑚) TVD (𝑚) (b) With regularization Figure 2.16: Inverted formation of model problem I using the inversion strategy of Example B.2, i.e., with input measurements corresponding to 65 logging positions per sample. 34 2 Solving Inverse Problems using Deep Learning 0 50 100 150 200 250 300 350 400 450 500 26 28 30 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (a) LWD coaxial measurement. Without regularization 0 50 100 150 200 250 300 350 400 450 500 26 28 30 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (b) LWD coaxial measurement. With regularization Figure 2.17: Model problem I. Comparison between F˝Iand F˝Iθ˚using the inversion strategy of Example B.2, i.e., with input measurements corresponding to 65 logging positions per sample. 35 2 Solving Inverse Problems using Deep Learning recovered model (see Figure 2.18b) is still quite different from the original one. Nonetheless, we observe a superb matching in the logs (see Figure 2.20), which indicates the presence of a different solution for the inverse problem. This confirms that the given measurements are insufficient to provide a unique solution for the inverse problem. For the case with regularization, inversion results (see Figure 2.18b) match the original model, and the corresponding logs properly approximate the synthetic ones, see Figure 2.21. Figures 2.22 and 2.23 confirm that our methodology delivers a proper training of the forward function approximation and the composition Fφ˚˝Iθ˚, respectively. 36 2 Solving Inverse Problems using Deep Learning 2.8.3.2 Model Problem II In this problem, we consider a 2.5m-thick conductive layer surrounded by two resistive layers. A well trajectory with a dip angle equal to 87˝crosses the formation. Figure 2.24 displays the original and predicted models by DL. This example illustrates some of the limitations of DNNs. In this case, the Earth models associated with part of the trajectory are outside the model problems considered in Section 2.1, which restrict to only one layer above and below the logging trajectory. Thus, the DNN is untrained for such models, and results cannot be trusted in those zones. Numerical results confirm these observations. Nonetheless, inaccurate inversion results are simple to identify by inspection of the logs (Figures 2.25 and 2.26). 0 20 40 60 80 100 120 140 160 180 200 220 240 40 42 44 46 HD (𝑚) TVD (𝑚) (a) Actual formation 0 20 40 60 80 100 120 140 160 180 200 220 240 40 42 44 46 HD (𝑚) TVD (𝑚) (b) Predicted formation using one logging position with regularization Figure 2.24: Model problem 2. Comparison between actual and predicted formations with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 43 2 Solving Inverse Problems using Deep Learning 0 20 40 60 80 100 120 140 160 180 200 220 240 26 28 30 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (a) LWD coaxial measurement 0 20 40 60 80 100 120 140 160 180 200 220 240 −49 −48.8 −48.6 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (b) Deep coaxial measurement 0 20 40 60 80 100 120 140 160 180 200 220 240 12 14 16 18 20 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (c) Geosignal measurement Figure 2.25: Model problem 2. Comparison between F˝Iand F˝Iθ˚with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 44 2 Solving Inverse Problems using Deep Learning 0 20 40 60 80 100 120 140 160 180 200 220 240 26 28 ◦vs 𝜙∗◦ HD (𝑚) Att. (𝑑𝐵) (a) LWD coaxial measurement 0 20 40 60 80 100 120 140 160 180 200 220 240 −49 −48.5 −48 ◦vs 𝜙∗◦ HD (𝑚) Att. (𝑑𝐵) (b) Deep coaxial measurement 0 20 40 60 80 100 120 140 160 180 200 220 240 10 15 20 ◦vs 𝜙∗◦ HD (𝑚) Att. (𝑑𝐵) (c) Geosignal measurement Figure 2.26: Model problem 2. Comparison between F˝Iand Fφ˚˝Iwith regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 45 2 Solving Inverse Problems using Deep Learning 2.8.3.3 Model Problem III We now consider a model formation exhibiting geological faults and two different well trajectories. For well trajectory 1, Figure 2.27 shows the model problem, logging trajectory, inversion results, and coaxial attenuation logs. Inversion results are excellent. When considering the second well trajectory shown in Figure 2.28, we observe good inversion results except at the proximity of points with horizontal distance (HD) equals to 75m and 350m. These inaccurate inversion results are easily identified by examination of the corresponding logs. 46 2 Solving Inverse Problems using Deep Learning 0 50 100 150 200 250 300 350 400 450 500 40 50 60 HD (𝑚) TVD (𝑚) (a) Actual formation 0 50 100 150 200 250 300 350 400 450 500 40 50 60 HD (𝑚) TVD (𝑚) (b) Predicted formation 0 50 100 150 200 250 300 350 400 450 500 26 28 30 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (c) LWD coaxial measurement Figure 2.27: Model problem III, trajectory 1. Comparison between actual and predicted formations and the corresponding coaxial logs with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 47 2 Solving Inverse Problems using Deep Learning 0 50 100 150 200 250 300 350 400 450 500 40 50 60 HD (𝑚) TVD (𝑚) (a) Actual formation 0 50 100 150 200 250 300 350 400 450 500 40 50 60 HD (𝑚) TVD (𝑚) (b) Predicted formation 0 50 100 150 200 250 300 350 400 450 500 26.5 27 27.5 ◦vs ◦𝜃∗ HD (𝑚) Att. (𝑑𝐵) (c) LWD coaxial measurement Figure 2.28: Model problem III, Trajectory 2. Comparison between actual and predicted formations and the corresponding coaxial logs with regularization using the inversion strategy of Example B.1, i.e., with input measurements corresponding to one logging position per sample. 48 3 Database Generation using IGA Deep Learning (DL) methods are fast, but require a massive training dataset. To decrease the online computational time during field operations, we often produce such a large dataset a priori (offline) using tens of thousands of simulations of borehole resistivity measurements (see [75]). To generate the database for DL inversion, we employ simulation methods to solve Maxwell’s equations with different conductivity distributions (Earth models). Since 3D simulations are expensive and possibly unaffordable when computing such large databases, it is common to reduce the Earth model dimensionality to two or one spatial dimensions using a Fourier or a Hankel transform. These transformations lead to the so-called 2.5D [3, 48, 92, 101, 135] and 1.5D [12, 100, 129] formulations, respectively. 1.5D simulations are inaccurate when dealing with geological faults. In this work, we focus on the efficient generation of a massive database using 2.5D simulations – as a preliminary stage for DL inversions. We propose the use of refined Isogeometric Analysis (rIGA) discretizations to generate databases for DL inversion of 2.5D geosteering electromagnetic (EM) measurements. 3.1 2.5D Variational Formulation of Electromagnetic Measurements 3.1.1 3D Wave Propagation Problem The two time-harmonic curl Maxwell’s equations describing the 3D wave propagation in an isotropic medium are ∇ˆE`iωµH“ ´iωµM ∇ˆH“ pσ`iωεqE(3.1) where Eis the electric field, His the magnetic field, iis the imaginary unit, σis the electric conductivity, εis the electric permittivity, µis the magnetic permeability, ω“2πf is the angular frequency, with fbeing the transmitter frequency, and Mis the time-harmonic magnetic source located at px0, y0, z0q and given by M“δpx´x0qδpy´y0qδpz´z0qrMx, My, Mzs P R3,(3.2) 49 3 Database Generation using IGA with δp¨q being the Dirac delta function defined as follows: δpx´x0q:“"8, x “x0, 0, x ‰x0.(3.3) To incorporate an integrable approximation of the Dirac delta function, we consider a bell-like representation for the delta function. For example, in the x direction, we approximate: δpx´x0q « 1 α?πexp „´´x´x0 α¯2,(3.4) where αis a positive value. From Maxwell’s equations, we obtain the following reduced wave formulation in terms of magnetic field H: $ ’ ’ & ’ ’ % Find H“ rHx, Hy, Hzs,with H:ΩĂR3ÑC3,such that: ∇ˆˆ1 σ`iωε∇ˆH˙`iωµH“ ´iωµM,in Ω, Eˆn“0,on BΩ, (3.5) where Ωis the domain of study and nis the unit normal (outward) vector on the boundary BΩ. We define Ωas a tensor-product box: Ω“ΩxˆΩyˆΩz“`´Lx{2, Lx{2˘ˆ`´Ly{2, Ly{2˘ˆ`´Lz{2, Lz{2˘, being Lx,Ly, and Lzpositive real constants. To introduce the weak formulation of this problem, we first define the Hpcurl; Ωqconforming functional spaces Hpcurl; Ωq:“ tW“ rWx, Wy, Wzs P pL2pΩqq3:∇ˆWP pL2pΩqq3u, H0pcurl; Ωq:“ tWPHpcurl; Ωq:Wˆn“0on BΩu.(3.6) The Hpcurl; Ωqspace is endowed with the inner product pW,HqHpcurl;Ωq:“ p∇ˆW,∇ˆHqpL2pΩqq3`pW,HqpL2pΩqq3 :“żΩp∇ˆWq˚¨p∇ˆHqdΩ`żΩ W˚¨HdΩ,(3.7) where ˚is the conjugate transpose of complex vector space and ¨denotes the inner product. 50 3 Database Generation using IGA We build the weak formulation by multiplying Eq. (3.5) with an arbitrary function WPH0pcurl; Ωq, using Green’s formula, and integrating over the domain Ω. The weak formulation is then $ & % Find HPH0pcurl; Ωq,such that for every WPH0pcurl; Ωq, ˆ∇ˆW,1 σ`iωε∇ˆH˙pL2pΩqq3`iωµpW,HqpL2pΩqq3“ ´iωµpW,MqpL2pΩqq3 (3.8) 3.1.2 2.5D Variational Formulation Herein, we focus on the case when the material properties are homogeneous along one spatial direction, e.g., y-axis. We denote the domain for this case as Ω:“ΩyˆΩx,z. We perform a Fourier transform along the y-axis to represent the 3D problem as a sequence of uncoupled 2D problems, one per Fourier mode. In this case, we define the magnetic field Has a series expansion using the complex exponentials: H:“ `8 ÿ β“´8 Hβexppi2πβy{Lyq,(3.9) where βis the Fourier mode number and Hβ““Hβ x, Hβ y, Hβ z‰with Hβ:Ωx,z Ă R2ÑC3. Fourier modes satisfy the following orthogonality relationships, where δi,j is the Kronecker delta: 1 LyżLy{2 ´Ly{2 exppi2πβ1y{Lyqexppi2πβ2y{Lyqdy “δβ1,β2.(3.10) By employing a test function of the form W:“1 Ly Wβexppi2πβy{Lyq,(3.11) and defining the Hpcurlβ;Ωx,zq-conforming functional spaces Hpcurlβ;Ωx,zq:“ Wβ“ rWβ x, Wβ y, Wβ zs P pL2pΩx,zqq3:Wβ yPH1pΩx,zq and ∇ˆrWβ x, Wβ zs P pL2pΩx,zqq2(, H0pcurlβ;Ωx,zq:“ WβPHpcurlβ;Ωx,zq:Wβˆn“0on BΩ(. (3.12) 51 3 Database Generation using IGA we build the following variational formulation from Eq. (3.5) by integrating over Ωx,z: $ ’ ’ & ’ ’ % Find H“ř`8 β“´8 Hβexppi2πβy{Lyq,HβPH0pcurl; Ωx,zq such that for every βPZand WβPH0pcurl; Ωx,zq, ˆ∇βˆWβ,1 σ`iωε∇βˆHβ˙pL2pΩx,zqq3`iωµpWβ,HβqpL2pΩx,z qq3“ ´iωµpWβ,MβqpL2pΩx,zqq3, (3.13) where ∇βˆWβ:“«iβ 2π Ly Wβ z´dWβ y dz ,dWβ x dz ´dWβ z dx ,dWβ y dx ´iβ 2π Ly Wβ xff,(3.14) and Mβis the time-harmonic magnetic source written in terms of the Fourier transform as Mβ“1 Ly δpx´x0qδpz´z0qrMx, My, Mzsexppi2πβy0{Lyq.(3.15) This formulation corresponds to the 2.5D variational formulation previously described by, e.g., [27] and [120]. Remark 3.1. To solve the variational problem of Eq. (3.13), we require an appropriate space in Ωx,z over which Wβ,∇ˆrWβ x,Wβ zs, and ∇Wβ yare integrable, i.e., rWβ x,Wβ zs P Hpcurl; Ωx,zqand Wβ yPH1pΩx,zq. Thus, we use the Hpcurlβ;Ωx,zq solution space – equivalent to the Hpcurl; Ωx,zq ˆ H1pΩx,zqmixed space – that fulfills the mentioned requirements. 3.2 Borehole Resistivity Measurement Acquisition System We consider the logging-while-drilling (LWD) instrument equipped with transmitters (Txi) and receivers (Rxj) of Figure 3.1. This tool is sensitive to resistivities within the range 0.2„500 Ω ¨m (phase resistivity) and 0.2„300 Ω ¨m (amplitude resistivity) under an operating frequency between 0.1 and 2 MHz [79]. For the sake of simplicity, herein, we restrict to two transmitters and two receivers symmetrically located around the tool center (see Figure 3.1) at an operating frequency of 2 MHz. Triaxial logging instruments generate measurements for all possible orientations of the transmitter–receiver pairs. We follow the notation presented by [34] 52 3 Database Generation using IGA z x l l h h L L Ωx,z T2R2R1T1Tool Trajectory Figure 3.4: A drawing of the computational domain Ωx,z and the tool trajectory. The central subgrid bounded by a magenta box is composed of a set of fine elements located in the proximity of the logging instrument. The remaining elements grow smoothly in size until reaching the boundary. Fourier modes. Due to the symmetry of the media along the ydirection, we only consider βP r0, Nfs. 3.5 Numerical Results In this section, we first assess the accuracy of the rIGA approach in a homogeneous medium. We also investigate the computational efficiency of the rIGA framework in comparison with IGA and FEM approaches. Then, we consider two model problems consisting of high-angle wells crossing spatially heterogeneous media with multiple geological faults. Finally, we produce our synthetic training dataset as a preliminary stage for DL inversion. In our simulations, we consider one operational mode of a commercial logging tool [158] with lR“10.16 cm and lT“56.8325 cm (see Figure3.1). We select the free space electric permittivity and magnetic permeability as ε“8.85ˆ10´12 F¨m´1and µ“4πˆ10´7N¨A´2, respectively. We also consider a transmitter frequency f“2 MHz. 3.5.1 Homogeneous Medium We assume the logging instrument is placed in a homogeneous medium with resistivity ρ“1{σ“100 Ω ¨m. This high-resistivity case is numerically more challenging than low-resistivity cases since it requires a larger number of Fourier 59 3 Database Generation using IGA modes and numerical precision. We consider a cube domain of length L“18 m for our first case study. 3.5.1.1 Accuracy Assessment To assess the accuracy and select certain discretization parameters, we compare the numerical attenuation ratio, A, and phase difference, P, given by Eq.(3.17) and Eq.(3.18), with the expected (i.e., exact) values, Aeand Pe, obtained from ρe“1{σ. Figure 3.5 shows the numerical errors, i.e., |1´A{Ae|and |1´P{Pe|, as a function of the number of Fourier modes when computing attenuation ratio and phase difference in a homogeneous medium. Herein, we select a domain with 64 ˆ64 elements to ensure a fast numerical solution for our measurements. We compare the results of the high-continuity Cp´1IGA with FEM and also with a rIGA discretization that employs 8 ˆ8 macroelements. This macroelement size provides the fastest results for moderate size domains (see [47]). We consider three different mesh sizes – h“0.025 m,0.033 m,and 0.050 m– and different polynomial degrees –p“3,4,and 5. The best results correspond to h“0.025 m (blue lines in the figure). We also observe that rIGA and FEM discretizations deliver lower errors compared to their IGA counterparts –when the same number of elements and polynomial degree are considered– taking into account that rIGA provides solutions with higher computational efficiency (see Section 3.5.1.2). To investigate the decay of the solution for each Fourier mode, we compare numerical results with the analytical 2.5D solution in the homogeneous medium presented by [120]. In particular, given Mzas the only nonzero component of the magnetic source, it is possible to analytically determine the coaxial magnetic field for each Fourier mode as follows: HZZ pβq“´iωµστβz`B2τβz Bz2,(3.34) with τβz“Mz 2π 1 Ly K0pCRqexppi2πβy0{Lyq,(3.35) where K0p¨q is the modified Bessel function of the second kind of order zero, and C“ p2πβ{Lyq2`iωµσ, (3.36) R“apx´x0q2`pz´z0q2.(3.37) 60 3 Database Generation using IGA 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m |1−A/Ae| IGA rIGA FEA 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m |1−A/Ae| IGA rIGA FEA 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m |1−A/Ae| IGA rIGA FEA 10 20 30 40 50 60 70 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m No. of Fourier modes, Nf |1−P/Pe| IGA rIGA FEA (a) p“3 10 20 30 40 50 60 70 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m No. of Fourier modes, Nf |1−P/Pe| IGA rIGA FEA (b) p“4 10 20 30 40 50 60 70 10−4 10−3 10−2 10−1 100 h=0.025 m h=0.033 m h=0.050 m No. of Fourier modes, Nf |1−P/Pe| IGA rIGA FEA (c) p“5 Figure 3.5: Numerical errors when computing attenuation ratio, A, and phase difference, P, in a homogeneous medium using IGA and rIGA discretizations, obtained by a 64ˆ64 element mesh with different element sizes hand polynomial degrees p. Figure 3.6 compares the decay of the numerical coaxial magnetic field HZZpβq with its analytical counterpart for some Fourier modes. Using a domain with 64 ˆ64 elements and h“0.025 m, we monitor the decay of the propagated waves at distances within the interval r0.2,1.0sm from the transmitters to ensure that the solutions at both receivers properly approximate the analytical ones. Results show that rIGA discretizations deliver increased accuracy for all tested polynomial degrees. 3.5.1.2 Computational Efficiency [46, 47] provide theoretical cost estimates of solving H1and Hpcurlqdiscrete spaces, respectively. Herein, we add these estimates to predict the cost of discretizing the Hpcurl; Ωx,zqˆH1pΩx,zqspace appearing in our 2.5D EM problem. We conclude that the cost of LU factorization of the rIGA matrix for this combined space is between Oppqand Opp2qtimes smaller than that for IGA. Details are omitted for the sake of simplicity. To numerically assess the computational efficiency confirming the aforementioned theoretical results, we consider two different grids in Ωx,z with 64 ˆ64 and 128 ˆ128 elements, respectively. Using continuity reduction, we split the 61 3 Database Generation using IGA 0.2 0.40.60.8 1 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 β=0 β=10 β=25 β=45 Distance from transmitter (m) |HZZ|(A ·m−1) Numerical (IGA) Numerical (rIGA) Analytical (a) p“3 0.2 0.40.60.8 1 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 β=0 β=10 β=25 β=45 Distance from transmitter (m) |HZZ|(A ·m−1) Numerical (IGA) Numerical (rIGA) Analytical (b) p“4 0.2 0.40.60.8 1 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 β=0 β=10 β=25 β=45 Distance from transmitter (m) |HZZ|(A ·m−1) Numerical (IGA) Numerical (rIGA) Analytical (c) p“5 Figure 3.6: Comparison of the decay of the numerical and analytical coaxial magnetic fields for some Fourier modes, obtained in a grid of 64 ˆ64 elements with h“0.025 m and different polynomial degrees. mesh symmetrically into macroelements whose sizes are powers of two. In this context, the maximum-continuity Cp´1IGA discretization is composed of one macroelement containing the entire grid, while C0FEM with minimum continuity across all element interfaces is composed of macroelements that contain only one element. Figure 3.7 shows the number of FLOPs and time required to solve the borehole resistivity problem for each Fourier mode per logging position. We compare the computational costs for different polynomial degrees and different continuity reduction levels of basis functions. The cost of rIGA reaches the minimum with 8 ˆ8 macroelements almost in all cases, confirming the theoretical estimates obtained from the results of [47]. Our numerical tests show that for a moderate size 2.5D problem, the reduction in the number of FLOPs is Oppqwith respect to IGA. When compared to FEM, rIGA delivers larger improvement factors. These improvement factors in terms of FLOPs also hold in terms of time when performing a sequential factorization. In our parallel PARDISO solver, we observe a small degradation of the rIGA improvement factors in terms of times in comparison to those obtained in terms of FLOPs (see Table 3.1). 62 3 Database Generation using IGA 1×1 2×2 4×4 8×8 16×16 32×32 64×64 1010 1011 1012 p=3 p=4 p=5 Cp−1IGA C0FEA Macroelement size FLOPs (a) FLOPs (64 ˆ64 elements) 1×1 2×2 4×4 8×8 16×16 32×32 64×64 10−1 100 101 p=3 p=4 p=5 Cp−1IGA C0FEA Macroelement size Time (s) (b) Time (64 ˆ64 elements) 1×1 2×2 4×4 8×8 16×16 32×32 64×64 128×128 1011 1012 1013 p=3 p=4 p=5 Cp−1IGA C0FEA Macroelement size FLOPs (c) FLOPs (128 ˆ128 elements) 1×1 2×2 4×4 8×8 16×16 32×32 64×64 128×128 100 101 102 p=3 p=4 p=5 Cp−1IGA C0FEA Macroelement size Time (s) (d) Time (128 ˆ128 elements) Figure 3.7: Computational cost in terms of FLOPs and time for solving a 2.5D borehole resistivity problem per logging position per Fourier mode. We test rIGA discretizations with two different grids of 64 ˆ64 and 128 ˆ128 elements. The computational times correspond to the use of parallel solver PARDISO using two threads. 63 3 Database Generation using IGA Table 3.1: Computational cost for the 2.5D borehole resistivity measurements per logging position per Fourier mode. We report the solution time and FLOPs when using Cp´1IGA, rIGA with 8 ˆ8 macroelements, and C0 FEM with the same number of elements and polynomial degree. The computational times correspond to the use of parallel solver PARDISO using two threads. Domain size Polynomial degree Discretization method Number of FLOPs Improvement factor (FLOPs) Time (s) Improvement factor (time) IGA 6.56e+10 0.407 3 rIGA 2.85e+10 IGA/rIGA 2.30 0.242 IGA/rIGA 1.68 FEA 1.65e+11 FEA/rIGA 5.79 0.999 FEA/rIGA 4.12 IGA 1.62e+11 0.875 64ˆ64 4 rIGA 4.97e+10 IGA/rIGA 3.26 0.419 IGA/rIGA 2.09 FEA 5.56e+11 FEA/rIGA 11.19 3.061 FEA/rIGA 7.29 IGA 3.10e+11 1.645 5 rGA 7.33e+10 IGA/rIGA 4.23 0.620 IGA/rIGA 2.65 FEA 1.37e+12 FEA/rIGA 18.69 6.806 FEA/rIGA 10.98 IGA 5.72e+11 3.144 3 rIGA 2.24e+11 IGA/rIGA 2.55 1.423 IGA/rIGA 2.21 FEA 1.31e+12 FEA/rIGA 5.85 7.456 FEA/rIGA 5.24 IGA 1.43e+12 6.903 128ˆ128 4 rIGA 3.56e+11 IGA/rIGA 4.02 2.495 IGA/rIGA 2.77 FEA 4.44e+12 FEA/rIGA 12.47 22.885 FEA/rIGA 9.17 IGA 2.82e+12 12.911 5 rGA 4.97e+11 IGA/rIGA 5.67 3.305 IGA/rIGA 3.91 FEA – FEA/rIGA – – FEA/rIGA – 64 3 Database Generation using IGA 3.5.2 Heterogeneous Media We further examine the accuracy of our rIGA approximation over two synthetic heterogeneous model problems. 3.5.2.1 One Geological Fault We consider the model problem of Figure 3.8 with a constant dip angle of 80˝. We consider the Logging While Drilling (LWD) instrument described in Section 3.2 and simulate measurements recorded over 200 equally-spaced logging positions throughout the well trajectory. ρ=50 Ω·m 0.0 0.5 1.0 1.5 2.0 0.0 7.5 15.0 ρ=3Ω·m ρ=1Ω·m Figure 3.8: Model problem with a constant dip angle of 80˝passing through a geological fault and three different materials (well trajectory is highlighted by a red dashed line). Dimensions are in meters. Figure 3.9 shows the apparent resistivities based on the attenuation ratio and phase difference (ρAand ρP, respectively). We obtain the results using Nf“70 and a rIGA discretization with 64ˆ64 elements, p“4, and 8ˆ8 macroelements. Results are in good agreement with those presented by [120]. 3.5.2.2 Two Geological Faults and Inclined Layers Figure 3.10 shows the second model problem containing two geological faults and inclined layers. The logging trajectory starts from a sandstone layer with a resistivity of ρ“3 Ω ¨m, and passes through an oil-saturated layer with ρ“100 Ω ¨m. The tool trajectory also passes through a water-saturated layer with ρ“0.5 Ω ¨m. In particular, inclined layers produce the so-called staircase approximations [25]. This phenomenon occurs because the physical interfaces of the conductivity model are not aligned with the element edges. Thus, the conductivity parameter takes different values inside some elements of the mesh. To tackle this issue, discretization techniques using nonfitting grids [26, 27] are available, but they have not been considered here for simplicity. 65 3 Database Generation using IGA 02468 10 12 14 16 100 101 102 ρA ρP ρe Horizontal distance (m) Resistivity (Ω·m) Figure 3.9: Apparent resistivities based on the attenuation ratio, ρA, and phase difference, ρP, for the first model problem, compared with the real (exact) resistivity, ρe. We obtain the results using a rIGA discretization with 64 ˆ64 elements, p“4, and 8 ˆ8 macroelements. 0.0 3.0 4.0 5.0 6.0 8.0 9.0 11.0 0.0 30.0 60.0 90.0 5◦ 5◦ 1◦ ρ=3Ω·m ρ=100 Ω·m ρ=0.5Ω·m Figure 3.10: Second model problem with two geological faults and inclined layers. The tool trajectory (red dashed line) has different dip angles and passes through sandstone (yellow), oil-saturated (gray), and watersaturated (green) layers. Dimensions are in meters. Figure 3.11 shows the apparent resistivities based on the attenuation ration and phase difference throughout the logging trajectory and compares their value with the exact resistivity. We simulate the resistivities at 1,080 logging positions with Nf“70. We use a rIGA discretization with 64 ˆ64 elements, p“4, and 8ˆ8 macroelements. 3.5.3 Database Generation for DL Inversion To produce our synthetic training dataset for DL inversion, we consider heterogeneous medium containing three different layers and six varying parameters at 66 3 Database Generation using IGA 0 10 20 30 40 50 60 70 80 90 100 101 102 ρA ρP ρe Horizontal distance (m) Resistivity (Ω·m) Figure 3.11: Apparent resistivities based on the attenuation ratio, ρA, and phase difference, ρP, for the second model problem, compared with the real (exact) resistivity, ρe. We employ a rIGA discretization with 64ˆ64 elements, p“4, and 8 ˆ8 macroelements. each logging position, as described in Figure 3.12 and Table 3.2. We select three different electrical conductivities: σcfor the central layer, and σuand σlfor the upper and lower layers, respectively. We assume the tool center is always within the middle layer and has vertical distances of duand dlfrom the upper and lower layers, respectively. The sixth varying parameter is the dip angle, ϕ, measured from the vertical direction. + Tool center du dl Tool trajectory ϕ σu σc σl Figure 3.12: Varying parameters at each logging position when producing the training dataset for DL inversion. In here, we create a dataset of 100,000 samples and compute the apparent resistivities obtained from random combinations within a given range of resistivities ρ“1{σP r1,100sΩ¨m (see Table 3.2). For generating the dataset, we use two different types of parallelization. One parallelization is related to the parallel 67 3 Database Generation using IGA Table 3.2: Varying parameters employed to generate the training dataset for DL inversion. Varying parameters Interval Electrical conductivity of the central layer log10pσc)r´2,0s Electrical conductivity of the upper layer log10pσu)r´2,0s Electrical conductivity of the lower layer log10pσl)r´2,0s Distance of the tool center from the upper layer log10pdu)r´2,1s Distance of the tool center from the lower layer log10pdl)r´2,1s Dip angle between the tool and the layered media ϕr80˝,100˝s factorization of the direct solver, and the other is the trivial parallelization based on scheduling the solutions of independent Earth models onto different processors. Using 40 threads, we solve for 20 different Earth models, each executing over two threads. Table 3.1 shows that the required time for matrix factorization of the 2.5D EM problem using optimal rIGA discretization with 64 ˆ64 grid, p“4, and 8 ˆ8 macroelements is about 0.42 seconds per Fourier mode. Considering Nf“70, and the additional time required for pre/postprocessing and inter-thread communications, each set of independent runs (consists of 20 different Earth models) takes about 40 seconds. Thus, we perform 5,000 sequential runs to construct our 100,000 samples in about 56 hours. To create a larger database, we could execute over a cluster of hundreds of CPUs/threads, expecting a perfect parallel scalability. Figure 3.13a depicts the graphs of attenuation ratio, A, versus phase difference, P, obtained from the 100,000 Earth models when using rIGA discretization for generating the database. Since there is a strong correlation between Aand P, the data distribution on the plot follows an almost straight line. We also display in Figure 3.13b the correlation between apparent resistivities based on attenuation ratio and phase difference. 68 4 Quadrature Rules when Solving PDEs using Deep Learning Simplifying terms we reach to ´żΩ vp∆u`fqdx “0,@vPV. (4.17) With the assumption that p∆u`fqis continuous that equality can only be satisfied if p∆u`fqpxq “ 0xPΩ,(4.18) and this means that uis solution of (4.1). To finish the proof, we show that the solution of (4.5) is unique. Suppose that u1PVand u2PVare solutions of (4.5), p∇u1,∇vq “ pf, vq`pg, vqΓN@vPV, p∇u2,∇vq “ pf, vq`pg, vqΓN@vPV. (4.19) Subtracting the equations from (4.19) and selecting v“u1´u2we obtain żΩp∇u1´∇u2q2dx “0,(4.20) which implies that ∇u1pxq´∇u2pxq “ ∇pu1´u2qpxq “ 0@xPΩ.(4.21) This implies that pu1´u2qpxqis constant on Ω, and with the boundary condition u1pxq “ 0 for xPΓDwe reach to u1pxq “ u2pxq,@xPΩ. 4.2.2 Least Squares Method Reordering the terms of Eq. (4.1a) and (4.1c), we define: "Gu:“u2`f x PΩ, Bu:“u1¨n´g x PΓN.(4.22) To introduce the Least Squares method, we define the function FLS :VÝÑ R, where the function vPVsatisfies the Dirichlet conditions: FLSpvq “ |pGv, Gvq|`|pBv, BvqΓN|.(4.23) We want to minimize the function FLSpvqsubject to the essential (Dirichlet) Boundary Conditions (BCs). We often find the minimum by taking the derivative equal to zero and ending up with a linear system of equations. In the context of 75 4 Quadrature Rules when Solving PDEs using Deep Learning DL, we can simply introduce the above loss function FLSpvqdirectly in our NN. Therefore, we want to find u“arg min vPV FLSpvq.(4.24) 4.3 Neural Network Implementation We train a NN, named uNN px;θq, with the following architecture. We define the trainable part of our NN with learnable parameters θ. We call it uθ. It is composed by: 1. An input layer. This layer receives the data in the form of a nˆdmatrix, where nis the number of samples and dis the dimension of the data. 2. One hidden dense layer with mneurons and a sigmoid activation function. 3. An output layer that delivers uθ. Then, we add non-trainable layers to our arquitecture in order to impose Equations (4.6) or (4.23). For that, we introduce: 4. A non-trainable layer to impose the Dirichlet boundary conditions. For that, we select a function φpxqthat satisfies the Dirichlet conditions of the problem and its value is nonzero everywhere else [76]. In this work, we select the following φpxqfunctions for 1D problems in the interval Ω “ ra, bs: φpxq “ ź xDPΓDpx´xDq.(4.25) Then, we generate a new output of the NN: uNN px;θq “ φpxquθpxqthat strongly imposes the homogeneous Dirichlet boundary conditions. 5. A non-trainable layer to compute the loss function FRor FLS following Eqs. (4.6) or (4.23). Within this layer, we evaluate the integrals and the derivatives. We consider different quadrature rules, being the quadrature points part of the input data of our NN, along with the physical points of the domain. For computing the derivatives, we use automatic differentiation, except in some specific cases, where we employ FDM. These cases are explicitly indicated throughout the text. 76 4 Quadrature Rules when Solving PDEs using Deep Learning Figure 4.1 shows a schematic graph of the described NN architecture. Our software is developed in Python and we use the library Tensorflow 2.0. To train the NN, we replace in Equations (4.7) or (4.24) the search space V by the manifold generated by our learnable parameters θincluded in our NN. The result of the minimization is a function uNN px, ˜ θq, where ˜ θare the optimal learnable parameters encountered as a result of the training. For simplicity, in the following we abuse notation and use the symbol uNN to denote also the solution uNN px, ˜ θqof our minimization problem. x. . .uθuNN F(uNN ) Impose BC Define loss Trainable hidden layer Figure 4.1: Sketch of the arquitecture of uNN . 4.4 Quadrature Rules We approximate our integrals from Eqs. (4.6) and (4.23) using a quadrature rule of the form żb a fpxqdx « n ÿ i“0 ωifpxiq,(4.26) where ωiare the weights and xiare the quadrature points. Examples of quadrature rules that follow the above formula include trapezoidal rule and Gaussian quadrature rules [90]. We classify these quadrature rules into two groups: (1) those that only employ points from the interior of the interval; and (2) those that evaluate the solution at one extreme point or more (a or b). Integration rules within the later group (e.g., the trapezoidal rule) are inadequate for our minimization problems because the integrand can be infinite at the boundary points in the case of singular solutions (e.g. model problem 1). Thus, we focus on 77 4 Quadrature Rules when Solving PDEs using Deep Learning quadrature rules that only evaluate the solution at interior points of the domain, with a special focus on Gaussian quadrature rules. 4.4.1 Illustration of Quadrature Problems in Neural Networks 4.4.1.1 Ritz Method We consider the two model problems from Section 4.1. We approximate upxq using the Ritz method. Thus, we search for a NN that minimizes the loss functional given by Eq. (4.6). Our NN has one hidden layer with 10 neurons (31 trainable weights). We use automatic differentiation to compute the derivatives and a three-point Gaussian quadrature rule to approximate the integrals within each element. We select the Stochastic Gradient Descent (SGD) optimizer. For model problem 1, we discretize our domain with four equal-size elements and execute 40,000 iterations during the optimization process. For model problem 2, we discretize our domain with ten equal-size elements and execute 200,000 iterations during the optimization process. Figures 4.2a and 4.2b describe the loss evolution of the training process. We obtain a lower loss than the optimum loss computed analytically using the exact solution (i.e., FRpuexactq). This has to be due to some numerical error, in this case, quadrature errors. 100101102103104 8 −1.54 −15 Epoch Loss FR(uexact) FR(uNN ) (a) Model problem 1. 100101102103104105 0 −666.67 −800 Epoch Loss FR(uexact) FR(uNN ) (b) Model problem 2. Figure 4.2: Loss evolution of the training process for our two model problems. Figures 4.3 and 4.4 compare the approximate and exact solutions. We observe a disastrous NN approximation due to quadrature errors. Figures 4.3b and 4.4b show that the gradient is (almost) zero at the training (quadrature) points of the first interval. Therefore, this value minimizes the numerical approximation of 78 4 Quadrature Rules when Solving PDEs using Deep Learning pu1, u1q. This behavior allows the approximated solution to reach larger values in the first interval, and consequently maximizing the term pf, uq, and minimizing the total loss: FRpuNN q “ 1 2pu1 NN , u1 NN q looooomooooon řqiωipu1 NN q2«0 ´pf, uNN q looomooon «8 ´pg, uNN qΓN» ´8 The described quadrature errors can be interpreted as overfitting over the derivative of the solution. 0 2 4 6 8 10 0 20 40 x u uNN uexact (a) Exact and approximate solutions. 012 20 30 40 50 Gauss points x u (b) Approximate solution in the interval r0,2.5s. The Gauss quadrature points corresponding to the first element are indicated in blue. Figure 4.3: Exact vs approximate Ritz method solutions of model problem 1 using four elements for evaluating FRpvqand a NN with 31 weights. 4.4.1.2 Least Squares Method We now consider the following one-dimensional problem: "´u2pxq “ 0xP p0,1q, up0q “ u1p1q “ 0,(4.27) where the exact solution is upxq “ 0. We can easily construct an approximating function uNN that satisfies Eq. (4.27) at the three considered Gaussian points and minimizes Eq. (4.23), while still being a poor approximation of the exact solution due to quadrature errors. Figure 4.5 shows an example. 79 4 Quadrature Rules when Solving PDEs using Deep Learning 0 2 4 6 8 10 0 50 100 150 x u uNN uexact (a) Exact and approximate solutions. 0 0.2 0.4 0.6 0.8 1 −20 0 20 40 Gauss points x u uNN uexact (b) Approximate solution in the interval r0,1s. The Gauss quadrature points corresponding to the first element are indicated in blue. Figure 4.4: Exact vs approximate Ritz method solutions of model problem 2 using ten elements for evaluation of FRpvqand a NN with 31 weights. 0 0.25 0.5 0.75 1 0 a uNN uexact Gauss points x u Figure 4.5: Exact (uexact “0) and approximated solution of a problem given by Eq. (4.27) and solved with the Least Square (LS) method. 80 4 Quadrature Rules when Solving PDEs using Deep Learning 4.5 Integral Approximation We now describe four different methods to improve the integral approximations. 4.5.1 Monte Carlo Integration We consider the following Monte Carlo integral approximation over a set of points xiP pa, bq, żb a fpxqdx «pb´aq n n ÿ i“1 fpxiq, xiP pa, bq @i“ t1,¨¨¨ , nu(4.28) In the above, points xiare randomly selected [1]. While this method is useful for high-dimensional integrals, for low dimensions (1D, 2D, 3D) the computational cost is high since the value of the integral approximation converges as 1{?n[150]. 4.5.2 Piecewise-polynomial Approximation We replace the original NN uNN by a piecewise-polynomial approximation u˚ NN,¨, where ¨represents the number of pieces of our piecewise-linear interpolator. This approximation can be exactly differentiated (e.q., via FDM) and integrated (via a Gaussian quadrature rule). Figure 4.6 shows an example when we train a NN and we build a piecewise-linear approximation of the NN with four elements. This method controls quadrature errors. However, it is inadequate for highdimensional problems as we need a mesh that is difficult to implement and integration becomes time consuming. 4.5.3 Adaptive Integration We first consider a training dataset over the interval pa, bqby taking an equidistant partition of nelements. Then, we define the validation set as a global h-refinement of the training dataset. Figure 4.7 shows an example of a training and the corresponding validation datasets. Then, for each element of the training mesh (e.g., E1in Figure 4.7), we compare the numerical integral over that element vs the sum of the integrals over the two corresponding elements on the validation dataset (in our case, E1 1`E2 1). If the integral values differ by more than a stipulated tolerance, we h-refine the training element, and we upgrade the validation dataset so it is built as a global h-refinement of the training dataset. This process is described in Algorithm 1. Figure 4.7 also shows a training and the corresponding validation dataset after refining the first and third elements. We are able to control the quadrature errors by adding new quadrature points to the training dataset. However, the simplest way to implement such method is 81 4 Quadrature Rules when Solving PDEs using Deep Learning 012345678910 0 2 4 6 x u uNN uexact u∗ NN,4 Figure 4.6: Neural Network approximation uNN and its piecewise-linear element approximation u˚ NN,4. Algorithm 1: Adaptive integration method Generate a training dataset; Generate the corresponding validation dataset; Set tolerance and maximum iteration number imax; while iăimax do for j“1,¨¨¨ , n do Compute integral values Ijover the training dataset of elements Ej; Compute integral values I1 j,I2 jover the validation dataset of elements E1 j,E2 j; if |pI1 j`I2 jq´Ij| ą then h-refine the Ej-th element of the training set; h-refine the E1 j-th and E2 j-th elements of the validation set; else continue; end end i“i`1; end 82 4 Quadrature Rules when Solving PDEs using Deep Learning h ••••• ab E1E2E3E4 (a) Training set h/2 ••••••••• ab E1 1E2 1E1 3E2 3 (b) Validation set Figure 4.7: Black points (dots) correspond to the original (a) training/(b) validation partitions and blue points (circles) are the points added by the refinement performed in the first and third elements. by using meshes, which posses a limitation on high-dimensional integrals. As an alternative to generating a mesh, one can randomly add points to the training set. This entails difficulties when designing an adaptive algorithm. In the same way that we propose an h-adaptive method, we can also work with p-adaptivity [14] or a combination of them (e.g., hp-adaptivity [38]). 4.5.4 Regularization Methods We now introduce a problem-specific regularizer designed to control the quadrature error. In a one-dimensional setting, we consider the integral functional FRas given in (4.6), and its approximation via a midpoint rule, ˆ FR, given by ˆ FRpuq “b´a N N ÿ j“1ˆ1 2|u1pxjq|2´fpxjqupxjq˙`gpbqupbq`gpaqupaq.(4.29) We note that g“0, except where the Neumann condition is imposed. While we focus on the Ritz method, a similar heuristic can be applied to the LS method. We introduce a function Rthat depends on the learnable parameters θof a given Neural Network uNN , such that for any Neural Network with a given architecture, |FRpuNN q´ ˆ FRpuNN q| ă Rpθq.(4.30) If we then consider a loss function Lgiven by Lpθq “ ˆ FRpuNN q`Rpθq,(4.31) 83 4 Quadrature Rules when Solving PDEs using Deep Learning then we may be able to improve the approximation of the quadrature rule, as the loss contains a term that by design controls the quadrature error. For simplicity, we consider only the case of a single-layer network with a onedimensional input, and the midpoint rule for calculating the integral over a uniform partition of pa, bq. We consider a mid-point rule as in (4.29), and define the interval length δ“b´a Nand intervals Ij“`δ 2`xj, xj`δ 2˘. We estimate the error of the midpoint rule to integrate Fas ˇˇˇˇˇżb a Fpxqdx ´ N ÿ j“1 Fpxjqδˇˇˇˇˇ“ˇˇˇˇˇ N ÿ i“1żxj`δ 2 xj´δ 2pFpxq´Fpxjqqdxˇˇˇˇˇ ď N ÿ i“1żxj`δ 2 xj´δ 2|Fpxq´Fpxjq|dx ď N ÿ i“1 max tPIj|F1ptq|żxj`δ 2 xj´δ 2|x´xj|dx “δ2 4 N ÿ i“1 max tPIj|F1ptq|. (4.32) This estimate scales as Op1 Nqfor fixed F, and thus, for a large number of integration points, we expect the estimate to be sufficiently accurate and to avoid “overdamping” of the loss. With (4.32) in mind, we estimate the local Lipschitz constants of the integrand as in (4.6). The numerical estimation of the Lipschitz constants of NNs has attracted attention, as they form a way of estimating the generalizability of a Neural Network, and have been used in the training process as a way to encourage accurate generalization [44, 51, 125]. As we are dealing with loss functions that involve derivatives of the Neural Network, we however need estimates of higher order derivatives of uNN . The approach that we employ is similar in spirit to the work of [88] for obtaining aposteriori error estimates in PINNs. Despite the arithmetic complications involved in calculating R, conceptually the idea reduces to an application of Taylor’s theorem. On a single interval of integration Ij, we have that for every x, there exists some ξxso that |F1pxq| “ |F1pxiq`px´xiqF2pξxq| ď |F1pxiq|` δ 2||F2||8.(4.33) We then find Rusing a combination of local and global estimates for the derivatives of the integrand corresponding to simple pointwise evaluations at the integration points and global estimates involving the Neural Network weights. The necessary steps are: 84 4 Quadrature Rules when Solving PDEs using Deep Learning 0 2 4 6 8 10 0 2 4 x u uNN uexact (a) Model problem 1. 0 2 4 6 8 10 0 50 100 x u uNN uexact (b) Model problem 2. Figure 4.11: Ritz method solution when using adaptive integration. 4.6.3 Regularization Methods We do not apply the regularization method to model problem 1 as the method requires sufficient regularity in order to provide the necessary estimates in the calculation of R. Since the solution is singular at x“0, the necessary Lipschitz bounds on the integral functional cannot be obtained within this framework. Instead, we aim to demonstrate that for problems that are sufficiently regular, our technique can avoid overfitting, and leave open the question as to how one may adapt the technique to singular problems for future work. We thus consider model problem 2. We propose the loss defined via Lpθq “ ˆ FRpuNN q`Rpθq.(4.51) Explicitly, ˆ FRpuNN q “ 10 N N ÿ j“1 1 2|u1 NN pxjq|2´2uNN pxjq´20uNN p20q,(4.52) where xj“10 N`i´1 2˘. 4.6.3.1 Experiment 1 We consider N“50 points, and a single layer network with M“10 neurons. We use the Adam optimizer with learning rate 10´2. We solve model problem 2 with two losses: with and without regularization. In both cases, we measure the metrics L,R, and ˆ FR. For validation, we use an equidistant partition of p0,10q 91 4 Quadrature Rules when Solving PDEs using Deep Learning with 49 points, so that we still use a midpoint rule but with different integration points. 0 2 4 6 8 10 0 50 100 x u uNN uexact (a) Exact and approximated solutions. 100101102103104 0 1 2 3 ·108 Epoch Loss Loss L (b) Evolution of Lduring training. 100101102103104 −600 −400 −200 0 Epoch Loss Loss validation Loss training (c) Evolution of ˆ FRduring training. 100101102103104 0 1 2 3 ·108 Epoch Loss Loss R (d) Evolution of Rduring training. Figure 4.12: The solution and training information for Experiment 1 without regularization. Figure 4.12 shows the results without regularization. As expected, we see in Figure 4.12a that the approximation is poor due to overfitting, which is most notable around x“0 and attained within 5000 epochs. Via the provided plots we can observe the beginning of overfitting in two distinct manners. First, we observe in Figure 4.12c that the value of ˆ FRevaluated over the validation data begins to diverge from the value on the training data, becoming apparent at around 1000 epochs. We also see this behaviour reflected in the evolution of Rin Figure 4.12d, with its most dramatic increase beginning around the same 92 4 Quadrature Rules when Solving PDEs using Deep Learning 0 2 4 6 8 10 0 50 100 x u uNN uexact (a) Exact and approximated solutions. 100101102103104 −600 −400 −200 0 Epoch Loss Loss validation Loss training (b) Evolution of Lduring training. 100101102103104 −600 −400 −200 0 Epoch Loss Loss ˆ FR (c) Evolution of ˆ FRduring training. 100101102103104 0 10 20 Epoch Loss Loss R (d) Evolution of Rduring training. Figure 4.13: The solution and training information for Experiment 1 with regularization. iteration. This rapid increase also provokes an increase in L, as seen in Figure 4.12b. This also indicates that even if Ris not used as part of the training process, its increase could be used as a metric to identify overfitting. Figure 4.13 describes the results with regularization and we observe a different behaviour. The approximation is generally good, and we do not see any signs of overfitting within 105epochs, as shown in Figure 4.13a. In particular, the values of ˆ FRat the training and validation data remain consistent in Figure 4.13c. Throughout Figure 4.13 we see that within 105epochs all metrics appear to have converged to a limiting value. We obtain final values L« ´644.22, ˆ FR« ´666.07, R«24.8. We recall that the true energy of the exact solution is FRpuexactq«´666.667, which suggests the quadrature rule is accurate. Notice that in the case without regularization, before overfitting became apparent, R 93 4 Quadrature Rules when Solving PDEs using Deep Learning had already attained values of around 1000, which is far larger than the value of Rat the obtained solution when regularization was used. 4.6.3.2 Experiment 2 We now consider a smaller N. As we expect Rto scale as 1 N, we anticipate a more adverse effect when Nis small. To view this, we consider the same problem of Experiment 1, where we now select N“20 integration points. We consider M“10 neurons and minimize our problem using the Adam optimizer with a learning rate of 10´2. As before, we consider the cases with and without regularization. 0 2 4 6 8 10 0 50 100 x u uNN uexact (a) Exact and approximated solutions. 100101102103 0 0.5 1 1.5 ·107 Epoch Loss Loss L (b) Evolution of Lduring training. 100101102103 −800 −600 −400 −200 0 Epoch Loss Loss validation Loss training (c) Evolution of ˆ FRduring training. 100101102103 0 0.5 1 1.5 ·107 Epoch Loss Loss R (d) Evolution of Rduring training. Figure 4.14: The solution and training information for Experiment 2 without regularization. 94 4 Quadrature Rules when Solving PDEs using Deep Learning 0 2 4 6 8 10 0 50 100 x u uNN uexact (a) Exact and approximated solutions. 100101102103104 −400 −200 0 Epoch Loss Loss validation Loss training (b) Evolution of Lduring training. 100101102103104 −600 −400 −200 0 Epoch Loss Loss ˆ FR (c) Evolution of ˆ FRduring training. 100101102103104 0 50 100 150 Epoch Loss Loss R (d) Evolution of Rduring training. Figure 4.15: The solution and training information for Experiment 2 with regularization. Figure 4.14 presents the loss evolution without regularization. We observe overfitting, which is accompanied by divergence of the loss on the validation dataset, as well as a rapid increase in R, with these features visible within 5000 epochs. Figure 4.15 presents the results with regularization. We observe no signs of overfitting, with the validation and training loss remaining close in Figure 4.15b. All metrics appear to have converged to a limiting value within 104epochs. However, the large value of Rat the found solution (approximately 140) has substantially changed the optimization problem so that the obtained minimizer is far from the desired solution. The final value of ˆ FRis around ´622, which is far from the desired value of ´666.67. This experiment highlights the fact that the regularizer becomes more effective when a large number of integration points are 95 4 Quadrature Rules when Solving PDEs using Deep Learning used. 96 5 Conclusions and Future Work 5.1 Conclusions In this dissertation, we first focus on the use of Deep Neural Networks (DNNs) for the inversion of borehole resistivity measurements for geosteering applications. We analyze the strong impact that different loss functions have on the prediction results. For this, we illustrate via a simple benchmark example that a traditional data misfit loss function delivers poor results. As a remedy, we propose the use of an Encoder-Decoder based or a two-step based loss function. These approaches generate two DNN approximations: one for the forward function and another one for the inverse operator. Then, we apply these two loss functions in a field example with synthetic data, and we obtain adequate results. To guarantee that the inverse DNN approximation provides meaningful results, we need to ensure that the training dataset contains sufficient samples. Otherwise, both forward and inverse DNN operators may provide incorrect solutions while still ensuring the composition of both operators is close to the identity. Thus, the approach is highly dependent on the existence of a sufficiently rich training dataset, which facilitates the learning process of the DNNs. To ensure that the inverse DNN approximation delivers significant results, we find it highly beneficial to add a regularization term to the loss function based on the existing training dataset. This reduces the richness we need to guarantee within the training datasets. Nevertheless, such regularization terms may hide alternative feasible solutions for the inverse operator, which may provide overconfidence in the results. Another possibility is to consider a two-step based loss function. Using this approach, we have shown that the inverse problem considered in this work admits different solutions that are physically feasible, a fact that was obscured when using the regularization term. Other critical limitations of DNNs we encounter in this work are: (a) the limited approximation capabilities of DNNs to reproduce discontinuous functions, (b) the need for a new dataset and trained DNN for each subsurface parametrization, and (c) the poor results they exhibit when they are evaluated over a sample that is outside the training dataset space. More importantly, it is often difficult to identify the source of poor results, which may include inadequate selections of: (i) loss function, (ii) DNN architecture, (iii) regularization term, (iv) train97 5 Conclusions and Future Work ing dataset, (v) optimization algorithm, (vi) rescaling operator and norms, (vii) model parameterization, (viii) approximation capabilities of DNNs, or simply (ix) the nature of the problem due to a lack of adequate measurements. To deal with the afforementioned sources of errors, we propose a careful step-by-step error control based on: (a) selecting adequate norms, (b) proper rescaling of the variables, (c) selecting a well suited loss function possibly with a regularization term, (d) analyzing the evolution of the different terms of the loss function, (e) studying multiple cross-plots of different nature, and (f) performing an in-depth assessment of the results over multiple realistic test examples. We also show it is possible to obtain a good-quality inversion of geosteering measurements with limited online computational cost, thus, suitable for real-time inversion. Moreover, the quality of the inversion results can be rapidly evaluated to detect its possible inaccuracies in the field and select alternative inversion methods when needed. As mentioned before, the DNN approximation of the inverse operator is highly dependent on the existence of a sufficiently rich training dataset. In the case of 1D layered formations, it is often feasible to produce the required dataset. However, for more complicated cases, for example, in 2D and 3D geometries, a direct extension may be limited due to the larger number of inversion variables and the extremely time-consuming process of producing an exhaustive dataset. Such a large database is essential for layer-by-layer estimation of the inverted Earth models, which may be used for real-time adjustments of the well trajectory during geosteering operations. In the second part of this dissertation, we propose the use of refined Isogeometric Analysis (rIGA) discretizations for generating a massive synthetic database for Deep Learning (DL) inversion of 2.5D borehole electromagnetic (EM) measurements. rIGA delivers computational savings of up to Oppqcompared to the high-continuity Isogeometric Analysis (IGA). When compared to a traditional Finite Element Method (FEM) with the same mesh size and polynomial degree, rIGA provides higher improvement factors. At the same time, rIGA provides sufficiently accurate solutions for geosteering purposes. To create a dataset for DL inversion, we first selected certain discretization parameters based on the results of several homogeneous solutions. Then, we checked the accuracy over homogeneous and heterogeneous media. Finally, we generated a synthetic database composed of 100,000 Earth models with the corresponding measurements in about 56 hours using a workstation equipped with two CPUs. In the last part of this work, we focus on the use of Neural Networks (NNs) for solving a Partial Differential Equation (PDE). We first illustrate via two simple examples how quadrature errors can destroy the quality of the 98 5 Conclusions and Future Work approximated solution when solving PDEs using DL methods. For this, we solve two simple 1D problems based on Poisson’s equation using the Deep Ritz Method (DRM) and a three-point Gaussian quadrature rule. Then, we propose four different alternatives to overcome the quadrature problems, discuss their advantages and limitations, and illustrate their performance. In high dimensions, Monte Carlo integration methods are the best choice. Regularizer methods are another option, but they are problem dependent and they need to be derived for each different architecture. Moreover, they require further analysis for highly nonlinear integrands. Furthermore, they are limited only to sufficiently smooth integral functionals. In addition, more complex NN architectures (which should be needed in higher dimensions) will hinder the derivation of R. In low dimensions (three or below), Monte Carlo integration is not competitive because of its low convergence speed. In these cases, adaptive integration exhibits faster convergence. In the cases of piecewise-linear approximation and regularizers, we are also able to overcome the quadrature problems, but the convergence speed is often slower and the accuracy is lower than with adaptive integration. 5.2 Future Work There are several possible future research lines regarding this work. The first one is to consider more complex Earth models, possibly containing geological faults or other relevant subsurface features, and analyze the performance of the Encoder-Decoder and two-step based loss functions. Another line of research consists of reducing the dataset size required for solving inverse borehole problems. Thus, decreasing the computational cost of creating the dataset. For this, we use using Active Learning techniques. Another option is to use Transfer Learning techniques for higher spatial dimensions, which can also alleviate data requirements to train the corresponding DNN. Concerning Chapter 4, one possible future work is to implement adaptive integration for 2D and 3D problems. In the same way, the piecewise-polynomial approximation could be improved by implementing r-adaptivity to optimize the grid. Ultimately, we aim at solving parametric PDEs using NNs. By this, we will be able to solve the forward problem corresponding to different Earth models using NNs. Thus, we will be able to properly and efficiently create the dataset needed to train our DNN that approximates the solution of inverse borehole problems. 99 6 Main Achievements 6.1 Scientific Achievements In the first part of this dissertation, we investigate appropriate loss functions to train a Deep Neural Network (DNN) when dealing with an inverse problem. In the second part of this work, we propose the use of refined Isogeometric Analysis (rIGA) discretizations to generate databases for DL inversion of 2.5D geosteering electromagnetic (EM) measurements. In the third part of this work, we analyze the problems associated with quadrature rules in Deep Learning (DL) methods when solving Partial Differential Equations (PDEs), and we propose several alternatives to overcome quadrature problems. 6.2 Peer-reviewed Publications 6.2.1 Journals 2022 J. A. Rivera, J. M. Taylor, ´ A. J. Omella and D. Pardo. On quadrature rules for solving Partial Differential Equations using Neural Networks. Computer Methods in Applied Mechanics and Engineering, 2022, vol. 393, p. 114710. https://doi.org/10.1016/j.cma.2022.114710 2021 A. Hashemian, D. Garcia, J. A. Rivera and D. Pardo. Massive database generation for 2.5 D borehole electromagnetic measurements using refined isogeometric analysis. Computers & Geosciences, 2021, vol. 155, p. 104808. https://doi.org/10.1016/j.cageo.2021.104808 2021 M. Shahriari, D. Pardo, J. A. Rivera, C. Torres-Verd´ın, A. Picon, J. Del Ser, S. Ossand´on and V. M. Calo. Error control and loss functions for the deep learning inversion of borehole resistivity measurements. International Journal for Numerical Methods in Engineering, 2021, vol. 122(6), p. 1629-1657. https://doi.org/10.1002/nme.6593 100 BIBLIOGRAPHY [40] R. Desbrandes and R. Clayton. Chapter 9 measurement while drilling. Developments in Petroleum Science, 38:251 – 279, 1994. (cited in page(s) 2) [41] C. Dupuis and J. M. Denichou. Automatic inversion of deep-directionalresistivity measurements for well placement and reservoir description. The Leading Edge, 34(5):504–512, 2015. (cited in page(s) 2) [42] W. Ee, J. Han, and A. Jentzen. Deep Learning-Based Numerical Methods for High-Dimensional Parabolic Partial Differential Equations and Backward Stochastic Differential Equations. To appear in Communications in Mathematics and Statistics, 5, 06 2017. (cited in page(s) 5) [43] A. Esteva, A. Robicquet, B. Ramsundar, V. Kuleshov, M. DePristo, K. Chou, C. Cui, G. Corrado, S. Thrun, and J. Dean. A guide to deep learning in healthcare. Nature Medicine, 25:24–29, 2019. (cited in page(s) 3) [44] M. Fazlyab, A. Robey, H. Hassani, M. Morari, and G. Pappas. Efficient and accurate estimation of Lipschitz constants for deep neural networks. Advances in Neural Information Processing Systems, 32:11427–11438, 2019. (cited in page(s) 84) [45] D. Fleisch. A student’s guide to Maxwell’s equations. Cambridge University Press, 2008. (cited in page(s) 4) [46] D. Garcia, D. Pardo, and V. M. Calo. Refined isogeometric analysis for fluid mechanics and electromagnetics. Computer Methods in Applied Mechanics and Engineering, 356:598–628, 2019. (cited in page(s) 5, 55, 56, 57, 61) [47] D. Garcia, D. Pardo, L. Dalcin, M. Paszy´nski, N. Collier, and V. M. Calo. The value of continuity: Refined isogeometric analysis and fast direct solvers. Computer Methods in Applied Mechanics and Engineering, 316:586–605, 2017. (cited in page(s) 4, 56, 57, 60, 61, 62) [48] S. Gernez, A. Bouchedda, E. Gloaguen, and D. Paradis. Aim4res, an opensource 2.5D finite differences MATLAB library for anisotropic electrical resistivity modeling. Computers & Geosciences, 135:104401, Feb. 2020. (cited in page(s) 4, 49) [49] M. Ghasemi, Y. Yang, E. Gildin, Y. Efendiev, and V. M. Calo. Fast multiscale reservoir simulations using pod-deim model reduction. Society of Petroleum Engineers, pages 1–18, 2015. (cited in page(s) 2) 107 BIBLIOGRAPHY [50] S. Goswami, C. Anitescu, S. Chakraborty, and T. Rabczuk. Transfer learning enhanced physics informed neural network for phase-field modeling of fracture. Theoretical and Applied Fracture Mechanics, 106:102447, 2020. (cited in page(s) 71) [51] H. Gouk, E. Frank, B. Pfahringer, and M. J. Cree. Regularisation of neural networks by enforcing lipschitz continuity. Machine Learning, 110(2):393– 416, 2021. (cited in page(s) 84) [52] A. G¨une¸s Baydin, B. A. Pearlmutter, A. Andreyevich Radul, and J. Mark Siskind. Automatic differentiation in machine learning: A survey. Journal of Machine Learning Research, 18:1–43, 2018. (cited in page(s) 6) [53] J. Gunning and M. E. Glinsky. Detection of reservoir quality using bayesian seismic inversion. Geophysics, 72(3):R37–R49, 2007. (cited in page(s) 3) [54] A. Gupta, A. Anpalagan, L. Guan, and A. S. Khwaja. Deep learning for object detection and scene perception in self-driving cars: Survey, challenges, and open issues. Array, 10:100057, 2021. (cited in page(s) 3) [55] J. Hadamard. Lectures on Cauchy’s problem in linear partial differential equations. Yale University Press, 1923. (cited in page(s) 7) [56] T. Hageman, K. M. P. Fathima, and R. de Borst. Isogeometric analysis of fracture propagation in saturated porous media due to a pressurised nonNewtonian fluid. Computers and Geotechnics, 112:272–283, Aug. 2019. (cited in page(s) 4) [57] J. Han, A. Jentzen, and W. Ee. Solving high-dimensional partial differential equations using deep learning. Proceedings of the National Academy of Sciences, 115, 07 2017. (cited in page(s) 5) [58] J. R. Hauser. Partial Differential Equations: The Finite Element Method. Numerical Methods for Nonlinear Engineering Models. Springer Netherlands, 2009. (cited in page(s) 4, 5) [59] K. He, X. Zhang, S. Ren, and J. Sun. Deep residual learning for image recognition. arXiv:1512.03385, 2015. (cited in page(s) 3, 16) [60] C. F. Higham and D. J. Higham. Deep learning: An introduction for applied mathematicians. Computing Research Repository, abs/1801.05894, 2018. (cited in page(s) 4, 16, 17) 108 BIBLIOGRAPHY [61] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194(39-41):4135– 4195, 2005. (cited in page(s) 4) [62] C. Hur´e, H. Pham, and X. Warin. Some machine learning schemes for highdimensional nonlinear PDEs. Math. Comput., 89:1547–1579, 2020. (cited in page(s) 5) [63] S. ichi Amari. Backpropagation and stochastic gradient descent method. Neurocomputing, 5(4):185–196, 1993. (cited in page(s) 6) [64] O. Ijasana, C. Torres-Verd´ın, and W. E. Preeg. Inversion-based petrophysical interpretation of logging-while-drilling nuclear and resistivity measurements. Geophysics, 78 (6):D473–D489, 2013. (cited in page(s) 2, 10) [65] A. D. Jagtap, E. Kharazmi, and G. E. Karniadakis. Conservative physicsinformed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 365:113028, 2020. (cited in page(s) 71) [66] B. Jan, H. Farman, M. Khan, M. Imran, I. U. Islam, A. Ahmad, S. Ali, and G. Jeon. Deep learning in big data analytics: A comparative study. Computers & Electrical Engineering, 75:275 – 287, 2019. (cited in page(s) 3) [67] K. H. Jin, M. T. McCann, E. Froustey, and M. Unser. Deep convolutional neural network for inverse problems in imaging. IEEE Transactions on Image Processing, 26(9):4509–4522, 2017. (cited in page(s) 16) [68] Y. Jin, X. Wu, J. Chen, Y. Huang, et al. Using a physics-driven deep neural network to solve inverse problems for lwd azimuthal resistivity measurements. In SPWLA 60th Annual Logging Symposium. Society of Petrophysicists and Well-Log Analysts, 2019. (cited in page(s) 19) [69] C. Johnson. Numerical Solution of Partial Differential Equations by the Finite Element Method. Dover Books on Mathematics Series. Dover Publications, Incorporated, 2012. (cited in page(s) 4, 73) [70] G. Karypis and V. Kumar. A fast and high quality multilevel scheme for partitioning irregular graphs. SIAM Journal on Scientific Computing, 20(1):359–392, 1998. (cited in page(s) 58) 109 BIBLIOGRAPHY [71] E. Kharazmi, Z. Zhang, and G. E. Karniadakis. VPINNs: Variational Physics-Informed Neural Networks For Solving Partial Differential Equations, 2019. (cited in page(s) 5, 6, 71) [72] E. Kharazmi, Z. Zhang, and G. E. Karniadakis. hp-vpinns: Variational physics-informed neural networks with domain decomposition. Computer Methods in Applied Mechanics and Engineering, 374:113547, 2021. (cited in page(s) 71) [73] R. Khodayi-Mehr and M. Zavlanos. Varnet: Variational neural networks for the solution of partial differential equations. In Proceedings of the 2nd Conference on Learning for Dynamics and Control, volume 120 of Proceedings of Machine Learning Research, pages 298–307, 2020. (cited in page(s) 71) [74] A. Krizhevsky, I. Sutskever, and G. E. Hinton. Imagenet classification with deep convolutional neural networks. In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 25, pages 1097–1105. Curran Associates, Inc., 2012. (cited in page(s) 16) [75] D. Kushnir, N. Velker, A. Bondarenko, G. Dyatlov, and Y. Dashevsky. Real-time simulation of deep azimuthal resistivity tool in 2D fault model using neural networks. In SPE Annual Caspian Technical Conference and Exhibition, Oct. 2018. (cited in page(s) 49) [76] I. Lagaris, A. Likas, and D. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. IEEE Transactions on Neural Networks, 9:987–1000, 1998. (cited in page(s) 76) [77] R. J. LeVeque. Finite difference methods for ordinary and partial differential equations: steady-state and time-dependent problems. SIAM, 2007. (cited in page(s) 4, 5) [78] D. B. Lindell, J. N. Martel, and G. Wetzstein. Autoint: Automatic integration for fast neural volume rendering. In Proceedings of the conference on Computer Vision and Pattern Recognition (CVPR), 2021. (cited in page(s) 6) [79] C. R. Liu. Theory of Electromagnetic Well Logging. Elsevier, Amsterdam, Netherlands, 2017. (cited in page(s) 52) 110 BIBLIOGRAPHY [80] L. O. Loseth and B. Ursin. Electromagnetic fields in planarly layered anisotropic media. Geophysical Journal International, 170:44–80, 2007. (cited in page(s) 10, 23) [81] L. Lu, X. Meng, Z. Mao, and G. Karniadakis. Deepxde: A deep learning library for solving differential equations. SIAM Review, 63:208–228, 2021. (cited in page(s) 5) [82] L. Lu, Y. Zheng, G. Carneiro, and L. Yang. Deep Learning for Computer Vision: Expert techniques to train advanced neural networks using TensorFlow and Keras. Springer, Switzerland, 2017. (cited in page(s) 3) [83] M. L¨angkvist, L. Karlsson, and A. Loutfi. A review of unsupervised feature learning and deep learning for time-series modeling. Pattern Recognition Letters, 42:11 – 24, 2014. (cited in page(s) 3) [84] Z. Ma, D. Liu, H. Li, and X. Gao. Numerical simulation of a multi-frequency resistivity logging-while-drilling tool using a highly accurate and adaptive higher-order finite element method. Advances in Applied Mathematics and Mechanics, 4(4):439–453, 2012. (cited in page(s) 4) [85] A. Malinverno and C. Torres-Verd´ın. Bayesian inversion of DC electrical measurements with uncertainties for reservoir monitoring. Inverse Problems, 16(5):1343–1356, oct 2000. (cited in page(s) 3) [86] Z. Mao, A. D. Jagtap, and G. E. Karniadakis. Physics-informed neural networks for high-speed flows. Computer Methods in Applied Mechanics and Engineering, 360:112789, 2020. (cited in page(s) 6, 71) [87] X. Meng, Z. Li, D. Zhang, and G. E. Karniadakis. Ppinn: Parareal physicsinformed neural network for time-dependent pdes. Computer Methods in Applied Mechanics and Engineering, 370:113250, 2020. (cited in page(s) 71) [88] S. Mishra and R. Molinaro. Estimates on the generalization error of physics informed neural networks (PINNs) for approximating PDEs. arXiv preprint arXiv:2006.16144, 2020. (cited in page(s) 5, 84) [89] D. Moghadas. One-dimensional deep learning inversion of electromagnetic induction data using convolutional neural network. Geophysical Journal International, 222(1):247–259, 2020. (cited in page(s) 16) [90] P. Moin. Fundamentals of Engineering Numerical Analysis. Cambridge University Press, 2010. (cited in page(s) 77) 111 BIBLIOGRAPHY [91] D. Mortari. Least-Squares Solution of Linear Differential Equations. Mathematics, 5, 02 2017. (cited in page(s) 73) [92] M. J. Nam, D. Pardo, and C. Torres-Verd´ın. Simulation of boreholeeccentered triaxial induction measurements using a Fourier hp finiteelement method. Geophysics, 78(1):D41–D52, Jan. 2013. (cited in page(s) 4, 49) [93] D. M. Nguyen, A. Evgrafov, and J. Gravesen. Isogeometric shape optimization for electromagnetic scattering problems. Progress In Electromagnetics Research B, 45:117–146, 2012. (cited in page(s) 4) [94] V. P. Nguyen, C. Anitescu, S. P. Bordas, and T. Rabczuk. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 117:89–116, 2015. (cited in page(s) 5) [95] C. M. B. Nunes and C. R´egis. GEMM3D: An edge finite element program for 3D modeling of electromagnetic fields and sensitivities for geophysical applications. Computers & Geosciences, 139:104477, June 2020. (cited in page(s) 4) [96] G. Pang, L. Lu, and G. E. Karniadakis. fPINNs: Fractional PhysicsInformed Neural Networks. SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 2019. (cited in page(s) 5) [97] D. Pardo, L. Demkowicz, C. Torres-Verd´ın, and M. Paszynski. Twodimensional high-accuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method. SIAM Journal on Applied Mathematics, 66(6):2085–2106, Jan. 2006. (cited in page(s) 4) [98] D. Pardo, P. J. Matuszyk, V. Puzyrev, C. Torres-Verdin, M. J. Nam, and V. M. Calo. Modeling of Resistivity and Acoustic Borehole Logging Measurements Using Finite Element Methods. Elsevier, 2021. (cited in page(s) 4) [99] D. Pardo and C. Torres-Verdin. Fast 1D inversion of logging-while-drilling resistivity measurements for the improved estimation of formation resistivity in high-angle and horizontal wells. Geophysics, 80 (2):E111–E124, 2014. (cited in page(s) 2, 4, 10) [100] D. Pardo and C. Torres-Verd´ın. Fast 1D inversion of logging-while-drilling resistivity measurements for improved estimation of formation resistivity in 112 BIBLIOGRAPHY high-angle and horizontal wells. Geophysics, 80(2):E111–E124, Mar. 2015. (cited in page(s) 49) [101] 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 the simulation of 3D alternating current borehole resistivity measurements. Computer Methods in Applied Mechanics and Engineering, 197(4548):3836–3849, Aug. 2008. (cited in page(s) 4, 49) [102] J. A. Parker, R. V. Kenyon, and D. E. Troxel. Comparison of interpolating methods for image resampling. IEEE Transactions on Medical Imaging, 2(1):31–39, 1983. (cited in page(s) 17) [103] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer. Automatic differentiation in pytorch. In NIPS-W, 2017. (cited in page(s) 3) [104] M. Paszy´nski, R. Grzeszczuk, D. Pardo, and L. Demkowicz. Deep learning driven self-adaptive hp finite element method. In Computational Science – ICCS 2021, pages 114–121. Springer International Publishing, 2021. (cited in page(s) 5) [105] C. G. Petra, O. Schenk, and M. Anitescu. Real-time stochastic optimization of complex energy systems on high-performance computers. Computing in Science & Engineering, 16(5):32–42, 2014. (cited in page(s) 58) [106] C. G. Petra, O. Schenk, M. Lubin, and K. G¨aertner. An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization. SIAM Journal on Scientific Computing, 36(2):C139–C162, 2014. (cited in page(s) 58) [107] L. Piegl and W. Tiller. The NURBS Book. Springer-Verlag, New York, NY, 2nd edition, 1997. (cited in page(s) 54) [108] S. Purushotham, C. Meng, Z. Che, and Y. Liu. Benchmarking deep learning models on large healthcare datasets. Journal of Biomedical Informatics, 83:112–134, 2018. (cited in page(s) 3) [109] V. Puzyrev. Deep learning electromagnetic inversion with convolutional neural networks. Geophysical Journal International, 218(2):817–832, 2019. (cited in page(s) 16) [110] T. M. Quan, D. G. C. Hildebrand, and W.-K. Jeong. Fusionnet: A deep fully residual convolutional neural network for image segmentation in connectomics, 2016. (cited in page(s) 16) 113 BIBLIOGRAPHY [111] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. A. Hamprecht, Y. Bengio, and A. Courville. On the Spectral Bias of Neural Networks, 2019. (cited in page(s) 5) [112] M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019. (cited in page(s) 6, 71) [113] M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics Informed Deep Learning (Part I): Data-driven Solutions of Nonlinear Partial Differential Equations. 2017. (cited in page(s) 5) [114] S. Ranjan and S. Senthamilarasu. Applied Deep Learning and Computer Vision for Self-Driving Cars: Build autonomous vehicles using deep neural networks and behavior-cloning techniques. Packt Publishing, 2020. (cited in page(s) 3) [115] J. N. Reddy. Introduction to the finite element method. McGraw-Hill Education, 2019. (cited in page(s) 4) [116] W. Ritz. ¨ Uber eine neue methode zur l¨osung gewisser variationsprobleme der mathematischen physik. Journal f¨ur die reine und angewandte Mathematik, 135:1–61, 1909. (cited in page(s) 73) [117] J. A. Rivera, D. Pardo, and E. Alberdi. Design of loss functions for solving inverse problems using deep learning. In Computational Science – ICCS 2020, pages 158–171. Springer International Publishing, 2020. (cited in page(s) 3) [118] J. A. Rivera, J. M. Taylor, ´ Angel J. Omella, and D. Pardo. On quadrature rules for solving partial differential equations using neural networks. Computer Methods in Applied Mechanics and Engineering, 393:114710, 2022. (cited in page(s) 5) [119] ´ A. Rodr´ıguez-Rozas and D. Pardo. A priori Fourier analysis for 2.5D finite elements simulations of logging-while-drilling (LWD) resistivity measurements. Procedia Computer Science, 80:782–791, 2016. (cited in page(s) 58) [120] ´ A. Rodr´ıguez-Rozas, D. Pardo, and C. Torres-Verd´ın. Fast 2.5D finite element simulations of borehole resistivity measurements. Computational Geosciences, 22(5):1271–1281, 2018. (cited in page(s) 4, 52, 53, 58, 60, 65) 114 BIBLIOGRAPHY [121] L. Ruthotto and E. Haber. Deep Neural Networks Motivated by Partial Differential Equations. Journal of Mathematical Imaging and Vision, 62, 2020. (cited in page(s) 5) [122] F. Sahli Costabal, Y. Yang, P. Perdikaris, D. E. Hurtado, and E. Kuhl. Physics-informed neural networks for cardiac activation mapping. Frontiers in Physics, 8, 2020. (cited in page(s) 71) [123] E. Samaniego, C. Anitescu, S. Goswami, V. Nguyen-Thanh, H. Guo, K. Hamdia, X. Zhuang, and T. Rabczuk. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 362:112790, 2020. (cited in page(s) 5) [124] A. Sarmiento, A. Cˆortes, D. Garcia, L. Dalcin, N. Collier, and V. Calo. PetIGA-MF: A multi-field high-performance toolbox for structurepreserving B-splines spaces. Journal of Computational Science, 18:117–131, 2017. (cited in page(s) 57) [125] K. Scaman and A. Virmaux. Lipschitz regularity of deep neural networks: analysis and efficient estimation. In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 3839–3848, 2018. (cited in page(s) 84) [126] O. Schenk and K. G¨artner. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems, 20(3):475–487, Apr. 2004. (cited in page(s) 58) [127] O. Schenk, K. G¨artner, and W. Fichtner. Efficient sparse LU factorization with left-right looking strategy on shared memory multiprocessors. BIT Numerical Mathematics, 40(1):158–176, 2000. (cited in page(s) 58) [128] D. J. Seifert, S. A. Dossary, R. E. Chemali, M. S. Bittar, A. A. Lotfy, J. L. Pitcher, and M. A. Bayrakdar. Deep electrical images, geosignal, and real-time inversion help guide steering decisions. Society of Petroleum Engineers, pages 1–9, 2009. (cited in page(s) 2) [129] M. Shahriari and D. Pardo. Borehole resistivity simulations of oil-water transition zones with a 1.5D numerical solver. Computational Geosciences, 24(3):1285–1299, Mar. 2020. (cited in page(s) 4, 49) 115 BIBLIOGRAPHY [130] M. Shahriari, D. Pardo, B. Moser, and F. Sobieczky. A deep neural network as surrogate model for forward simulation of borehole resistivity measurements. Procedia Manufacturing, 42:235 – 238, 2020. International Conference on Industry 4.0 and Smart Manufacturing (ISM 2019). (cited in page(s) 4, 13) [131] M. Shahriari, D. Pardo, A. Pic´on, A. Galdran, J. D. Ser, and C. TorresVerd´ın. A deep learning approach to the inversion of borehole resistivity measurements. Computational Geosciences, 2020. (cited in page(s) 2, 4, 13) [132] M. Shahriari, D. Pardo, J. Rivera, C. Torres-Verd´ın, A. Picon, J. Del Ser, S. Ossand´on, and V. Calo. Error control and loss functions for the deep learning inversion of borehole resistivity measurements. International Journal for Numerical Methods in Engineering, 122, 2020. (cited in page(s) 3) [133] M. Shahriari, S. Rojas, D. Pardo, A. Rodr´ıguez-Rozas, S. A. Bakr, V. M. Calo, and I. Muga. A numerical 1.5D method for the rapid simulation of geophysical resistivity measurements. Geosciences, 8(6):1–28, 2018. (cited in page(s) 2, 4, 8, 10) [134] S. Shahrokhabadi, T. D. Cao, and F. Vahedifard. Isogeometric analysis through B´ezier extraction for thermo-hydro-mechanical modeling of saturated porous media. Computers and Geotechnics, 107:176–188, Mar. 2019. (cited in page(s) 4) [135] J. Shen and W. Sun. 2.5-D modeling of cross-hole electromagnetic measurement by finite element method. Petroleum Science, 5(2):126–134, May 2008. (cited in page(s) 4, 49) [136] K. Shukla, P. C. Di Leoni, J. Blackshire, D. Sparkman, and G. E. Karniadakis. Physics-informed neural network for ultrasound nondestructive quantification of surface breaking cracks. Journal of Nondestructive Evaluation, 39(3):61, 2020. (cited in page(s) 71) [137] A. Simona, L. Bonaventura, C. de Falco, and S. Sch¨ops. IsoGeometric approximations for electromagnetic problems in axisymmetric domains. Computer Methods in Applied Mechanics and Engineering, 369:113211, Sept. 2020. (cited in page(s) 4) [138] R. N. Simpson, Z. Liu, R. V´azquez, and J. A. Evans. An isogeometric boundary element method for electromagnetic scattering with compatible B-spline discretizations. Journal of Computational Physics, 362:264–289, June 2018. (cited in page(s) 4) 116