Numerical solution of 3-D electromagnetic problems in exploration geophysics and its implementation on massively parallel computers

The growing significance, technical development and employment of electromagnetic (EM) methods in exploration geophysics have led to the increasing need for reliable and fast techniques of interpretation of 3-D EM data sets acquired in complex geological environments. The first and most important step to creating an inversion method is the development of a solver for the forward problem. In order to create an efficient, reliable and practical 3-D EM inversion, it is necessary to have a 3-D EM modelling code that is highly accurate, robust and very fast. This thesis focuses precisely on this crucial and very demanding step to building a 3-D EM interpretation method. The thesis presents as its main contribution a highly accurate, robust, very fast and extremely scalable numerical method for 3-D EM modelling in geophysics that is based on finite elements (FE) and designed to run on massively parallel computing platforms. Thanks to the fact that the FE approach supports completely unstructured tetrahedral meshes as well as local mesh refinements, the presented solver is able to represent complex geometries of subsurface structures very precisely and thus improve the solution accuracy and avoid misleading artefacts in images. Consequently, it can be successfully used in geological environments of arbitrary geometrical complexities. The parallel implementation of the method, which is based on the domain decomposition and a hybrid MPI-OpenMP scheme, has proved to be highly scalable - the achieved speed-up is close to the linear for more than a thousand processors. Thanks to this, the code is able to deal with extremely large problems, which may have hundreds of millions of degrees of freedom, in a very efficient way. The importance of having this forward-problem solver lies in the fact that it is now possible to create a 3-D EM inversion that can deal with data obtained in extremely complex geological environments in a way that is realistic for practical use in industry. So far, such imaging tool has not been proposed due to a lack of efficient, parallel FE solutions as well as the limitations of efficient solvers based on finite differences. In addition, the thesis discusses physical, mathematical and numerical aspects and challenges of 3-D EM modelling, which have been studied during my research in order to properly design the presented software for EM field simulations on 3-D areas of the Earth. Through this work, a physical problem formulation based on the secondary Coulomb-gauged EM potentials has been validated, proving that it can be successfully used with the standard nodal FE method to give highly accurate numerical solutions. Also, this work has shown that Krylov subspace iterative methods are the best solution for solving linear systems that arise after FE discretisation of the problem under consideration. More precisely, it has been discovered empirically that the best iterative method for this kind of problems is biconjugate gradient stabilised with an elaborate preconditioner. Since most commonly used preconditioners proved to be either unable to improve the convergence of the implemented solvers to the desired extent, or impractical in the parallel context, I have proposed a preconditioning technique for Krylov methods that is based on algebraic multigrid. Tests for various problems with different conductivity structures and characteristics have shown that the new preconditioner greatly improves the convergence of different Krylov subspace methods, which significantly reduces the total execution time of the program and improves the solution quality. Furthermore, the preconditioner is very practical for parallel implementation. Finally, it has been concluded that there are not any restrictions in employing classical parallel programming models, MPI and OpenMP, for parallelisation of the presented FE solver. Moreover, they have proved to be enough to provide an excellent scalability for it

Computationally Efficient Steady--State Simulation Algorithms for Finite-Element Models of Electric Machines.

The finite element method is a powerful tool for analyzing the magnetic characteristics of electric machines, taking account of both complex geometry and nonlinear material properties. When efficiency is the main quantity of interest, loss calculations can be affected significantly due to the development of eddy currents as a result of Faraday’s law. These effects are captured by the periodic steady-state solution of the magnetic diffusion equation. A typical strategy for calculating this solution is to analyze an initial value problem over a time window of sufficient length so that the transient part of the solution becomes negligible. Unfortunately, because the time constants of electric machines are much smaller than their excitation period at peak power, the transient analysis strategy requires simulating the device over many periods to obtain an accurate steady-state solution. Two other categories of algorithms exist for directly calculating the steady-state solution of the magnetic diffusion equation; shooting methods and the harmonic balance method. Shooting methods search for the steady-state solution by solving a periodic boundary value problem. These methods have only been investigated using first order numerical integration techniques. The harmonic balance method is a Fourier spectral method applied in the time dimension. The standard iterative procedures used for the harmonic balance method do not work well for electric machine simulations due to the rotational motion of the rotor. This dissertation proposes several modifications of these steady-state algorithms which improve their overall performance. First, we demonstrate how shooting methods may be implemented efficiently using Runge-Kutta numerical integration methods with mild coefficient restrictions. Second, we develop a preconditioning strategy for the harmonic balance equations which is robust against large time constants, strong nonlinearities, and rotational motion. Third, we present an adaptive framework for refining the solutions based on a local error criterion which further reduces simulation time. Finally, we compare the performance of the algorithms on a practical model problem. This comparison demonstrates the superiority of the improved steady-state analysis methods, and the harmonic balance method in particular, over transient analysis.PhDElectrical Engineering: SystemsUniversity of Michigan, Horace H. Rackham School of Graduate Studieshttp://deepblue.lib.umich.edu/bitstream/2027.42/113322/1/pries_1.pd

Overview of Large-Scale Computing: The Past, the Present, and the Future

published_or_final_versio

Three-dimensional modelling and inversion of controlled source electromagnetic data

The marine Controlled Source Electromagnetic (CSEM) method is an important and almost self-contained discipline in the toolkit of methods used by geophysicists for probing the earth. It has increasingly attracted attention from industry during the past decade due to its potential in detecting valuable natural resources such as oil and gas. A method for three-dimensional CSEM modelling in the frequency domain is presented. The electric field is decomposed in primary and secondary components, as this leads to a more stable solution near the source position. The primary field is computed using a resistivity model for which a closed form of solution exists, for example a homogeneous or layered resistivity model. The secondary electric field is computed by discretizing a second order partial differential equation for the electric field, also referred in the literature as the vector Helmholtz equation, using the edge finite element method. A range of methods for the solution of the linear system derived from the edge finite element discretization are investigated. The magnetic field is computed subsequently, from the solution for the electric field, using a local finite difference approximation of Faraday’s law and an interpolation method. Tests, that compare the solution obtained using the presented method with the solution computed using alternative codes for 1D and 3D synthetic models, show that the implemented approach is suitable for CSEM forward modelling and is an alternative to existing codes. An algorithm for 3D inversion of CSEM data in the frequency domain was developed and implemented. The inverse problem is solved using the L-BFGS method and is regularized with a smoothing constraint. The inversion algorithm uses the presented forward modelling scheme for the computation of the field responses and the adjoint field for the computation of the gradient of the misfit function. The presented algorithm was tested for a synthetic example, showing that it is capable of reconstructing a resistivity model which fits the synthetic data and is close to the original resistivity model in the least-squares sense. Inversion of CSEM data is known to lead to images with low spatial resolution. It is well known that integration with complementary data sets mitigates this problem. It is presented an algorithm for the integration of an acoustic velocity model, which is known a priori, in the inversion scheme. The algorithm was tested in a synthetic example and the results demonstrate that the presented methodology is promising for the improvement of resistivity models obtained from CSEM data

An iterative semi-implicit scheme with robust damping

An efficient, iterative semi-implicit (SI) numerical method for the time integration of stiff wave systems is presented. Physics-based assumptions are used to derive a convergent iterative formulation of the SI scheme which enables the monitoring and control of the error introduced by the SI operator. This iteration essentially turns a semi-implicit method into a fully implicit method. Accuracy, rather than stability, determines the timestep. The scheme is second-order accurate and shown to be equivalent to a simple preconditioning method. We show how the diffusion operators can be handled so as to yield the property of robust damping, i.e., dissipating the solution at all values of the parameter \mathcal D\dt, where $\mathcal D$ is a diffusion operator and \dt the timestep. The overall scheme remains second-order accurate even if the advection and diffusion operators do not commute. In the limit of no physical dissipation, and for a linear test wave problem, the method is shown to be symplectic. The method is tested on the problem of Kinetic Alfv\'en wave mediated magnetic reconnection. A Fourier (pseudo-spectral) representation is used. A 2-field gyrofluid model is used and an efficacious k-space SI operator for this problem is demonstrated. CPU speed-up factors over a CFL-limited explicit algorithm ranging from $\sim20$ to several hundreds are obtained, while accurately capturing the results of an explicit integration. Possible extension of these results to a real-space (grid) discretization is discussed.Comment: Submitted to the Journal of Computational Physics. Clarifications and caveats in response to referees, numerical demonstration of convergence rate, generalized symplectic proo

Fast iterative boundary element methods for high-frequency scattering problems in 3D elastodynamics

International audienceThe fast multipole method is an efficient technique to accelerate the solution of large scale 3D scattering problems with boundary integral equations. However, the fast multipole accelerated boundary element method (FM-BEM) is intrinsically based on an iterative solver. It has been shown that the number of iterations can significantly hinder the overall efficiency of the FM-BEM. The derivation of robust preconditioners for FM-BEM is now inevitable to increase the size of the problems that can be considered. The main constraint in the context of the FM-BEM is that the complete system is not assembled to reduce computational times and memory requirements. Analytic preconditioners offer a very interesting strategy by improving the spectral properties of the boundary integral equations ahead from the discretization. The main contribution of this paper is to combine an approximate adjoint Dirichlet to Neumann (DtN) map as an analytic preconditioner with a FM-BEM solver to treat Dirichlet exterior scattering problems in 3D elasticity. The approximations of the adjoint DtN map are derived using tools proposed in [40]. The resulting boundary integral equations are preconditioned Combined Field Integral Equations (CFIEs). We provide various numerical illustrations of the efficiency of the method for different smooth and non smooth geometries. In particular, the number of iterations is shown to be completely independent of the number of degrees of freedom and of the frequency for convex obstacles