Fast Solvers for Time-Harmonic Maxwell's Equations in 3D by Dhavide Arjunan Aruliah BSc. Simon Eraser University 1993 MSc. Simon Eraser University 1996 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E O F D o c t o r of Ph i losophy in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Computer Science) We accept this thesis as conforming to the required standard The University of British Columbia August 2001 (c) Dhavide Arjunan Aruliah, 2001 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstrac t The speed of iterative solvers for discretizations of partial differential equations (PDEs) is a significant bottleneck in the performance of codes de-signed to solve large-scale electromagnetic inverse problems. A single data in-version requires solving Maxwell's equations dozens if not hundreds of times. An inherent difficulty in geophysical contexts is that the conductivity and permeability coefficients may exhibit discontinuities spanning several orders of magnitude. Furthermore, in the air, the conductivity effectively vanishes. In standard formulations of Maxwell's equations, the curl operator that dom-inates the P D E operator leads to strong mixing of field components and ill-conditioning of linear systems resulting from standard discretizations. The primary objective of this research is to build fast iterative solvers for the forward-modeling problem associated with electromagnetic inverse prob-lems in the frequency domain. Toward this goal, a Helmholtz decomposition of the electric field using a Coulomb gauge condition recasts the P D E problem in terms of scalar and vector potentials. The resulting indefinite system is then stabilized by addition of a vanishing term that lies in the kernel of the domi-nant curl operator. Finally, an extra differentiation recasts the P D E system in a diagonally-dominant form reminiscent of a "pressure-Poisson" formulation for incompressible fluid flow. The continuous P D E problem obtained is equiv-alent to the original Maxwell's system but has a structure that is amenable to reliable solution techniques. Using a finite-volume scheme, the P D E is discretized on a staggered grid in three dimensions. The discretization obtained possesses conservation properties typical of finite-volume methods. Furthermore, interface conditions imposed by discontinuities in the material coefficients are sensibly accounted for in deriving the discretization. Although the simple representation of the ii media on a Cartesian tensor-product grid uses staircase approximations of surfaces of discontinuity of the material coefficients, some analysis and a nu-merical study demonstrate the suitability of such coarse approximations for diffusive problems. The discretization yields a non-Hermitian sparse linear system of al-gebraic equations; various preconditioners for Krylov-subspace methods are described, analyzed, implemented, and tested. Of particular interest is a multigrid preconditioner that exploits both the structure of the P D E prob-lem and the availability of well-established solvers for elliptic P D E problems (in particular, Dendy's B O X M G solver). The end result is a robust solver for the forward-modeling equations that can be incorporated within a competitive inverse problem code. iii Contents Abstract ii Contents iv List of Tables vii List of Figures viii Acknowledgements ix Dedication x 1 Introduction 1 1.1 Background Sketch 3 1.2 Outline of research 9 2 Background Theory of Electromagnetism 14 2.1 Maxwell's Equations 14 2.2 Derived Models 16 2.2.1 Constitutive Relations 17 iv 2.2.2 Stationary Models: Electrostatics and Magnetostatics . 19 2.2.3 The Time-Harmonic Model 21 2.2.4 The Quasistatic Model 22 2.2.5 Electric Source Currents 23 2.2.6 Magnetic Source Currents 25 2.2.7 The Magnetotelluric Model 26 2.3 Boundary Conditions 29 2.3.1 Frequency-Domain Boundary Conditions 31 2.3.2 Interface Conditions 33 2.4 Second-order P D E Formulations 36 2.5 Weak Formulations 38 3 Vector and Scalar Potential Formulations 43 3.1 The Forward-modeling Problem 44 3.1.1 Helmholtz Decomposition 47 3.1.2 Stabilization 49 3.1.3 "Pressure-Poisson" Formulation 50 3.2 Boundary Conditions 53 3.2.1 Interface Conditions 59 3.3 Weak Formulations 61 4 Finite-Volume Discretizations 63 4.1 Domain of Discretization 63 4.2 The Yee Discretization 66 4.2.1 The Discrete Fields 68 v 4.2.2 Derivation of Yee Scheme 70 4.3 Discretization Using Potentials 74 4.3.1 The Discrete Fields 76 4.3.2 Derivation of the Finite-Volume Scheme 78 4.3.3 The Discrete Linear System 84 4.4 Grid Effects 87 4.4.1 Model Problem Analysis 89 4.4.2 Numerical Study 95 5 Construction and Analysis of Iterative Solvers 103 5.1 Iterative Methods 104 5.1.1 Preconditioned Krylov-Subspace Methods 104 5.2 Analytical Framework 110 5.2.1 Discrete Operators on a Periodic Grid I l l 5.3 Spectral Bounds 114 5.4 Numerical Experiments 121 6 Conclusions 132 6.1 Research Summary 132 6.2 Future Research 136 Bibliography 139 vi Lis t of Tables 2.1 Field quantities in Maxwell's equations 15 2.2 Constitutive parameters for Maxwell's equations 18 4.1 a = a-ns- results with vertical shifts 98 4.2 a = <JE'- results with vertical shifts 100 4.3 cr = <7E: results with fine & coarse resolutions 102 5.1 cr = a-Q and electric dipole source on uniform grids 125 5.2 a = <7B a n d electric dipole source on nonuniform grids 126 5.3 cr = O B and magnetic dipole source on uniform grids 127 5.4 a = CTB and magnetic dipole source on nonuniform grids. . . . 127 5.5 a = <TE and electric dipole source on uniform grids 128 5.6 cr = O B over a range of frequencies 128 vii Lis t of Figures 2.1 Infinitesimal volume V§ for interface conditions 34 3.1 A typical geophysical scenario 45 4.1 Primary and dual grids for discretization 66 4.2 Staggering of components of Eh for Yee scheme 70 4.3 Staggering of components of Hh for Yee scheme 71 4.4 Cross-sections of a, a, and the support Sh of a — a 91 4.5 Cross-sections of <THS and <7E 97 4.6 Cross-sections of <7E resolved on fine and coarse grids 99 5.1 Eigenvalue range and condition number vs. x H9 5.2 Cross-section of <XB 124 5.3 Comparison of A M G and MM 131 viii Acknowledgements I have many people to thank. To start, the staff and members of the Dept. of Computer Science and the Institute of Applied Mathematics at U B C have been wonderfully helpful. I offer thanks also to my friends; too many to list, they know who they are. And of course, my family has always given me unconditional love and support without which I could not have gotten so far. My graduate experiences have taught me much about scientific inquiry and collaboration. I want to thank the members of the Scientific Computation and Visualization group as well as the UBC-GIF for greatly broadening my perspective. Drs. David Moulton and Joel Dendy assisted me greatly with the generous loan of their codes, as did Dr. Yair Shapira. I thank Drs. Doug Oldenburg, Jim Varah, and Brian Wetton for their helpful advice as members of my supervisory committee. I am also deeply grateful for the unbounded enthusiasm and scientific insight of my friend and colleague Dr. Eldad Haber. It is encouraging that the spirit of cooperation and camaraderie thrives in my scientific peers in spite of other pressures. Finally, I am extremely indebted to my research supervisor, Dr. Uri Ascher. He has set an example as a scientific researcher, as a teacher, and as a person that will be difficult to live up to in my future career. His support, guidance, inspiration, and friendship have been essential through difficult times over the last few years. Thank you, Uri, for everything. D H A V I D E A R J U N A N A R U L I A H The University of British Columbia August 2001 ix To my mother's parents, Daniel Chellathurai Arulanantham and Grace Emily Athisayam Arulanantham, and my father's parents, Vethakutty Charles Aruliah and Clara Arunothayam Sinnathangam Aruliah. I wish I had been able to know you better and I hope have made you proud. x Chapter 1 I n t r o d u c t i o n A crucial bottleneck in the performance of computational inverse prob-lems is the efficiency of available solvers for the associated forward-modeling problems. Generic strategies for parameter function estimation are based on an iteration within which the forward-modeling problem is solved using the current parameters, the obtained solution is somehow used to update the pa-rameters, and the process is repeated until the parameter function converges in some suitable norm [40]. For large-scale inverse problems, inversion of the forward-modeling problem within each iteration of the inverse problem gener-ally requires some iterative algorithm in itself. Since, in principle, dozens if not hundreds of forward-modeling solves are needed for a single data inversion, the embedded forward-modelling solver must be quick and robust. Within the geophysical community, an important family of inverse prob-lems used for prospecting derives from electromagnetic surveys [2, 42, 52, 75, 83, 85, 86]. In particular, the electromagnetic response due to some ground-based or aerial transmitter is measured over some terrain. The basic geophysi-1 cal inverse problem considered, then, is to construct a model of the conductiv-ity profile under the terrain based on the knowledge of the applied electromag-netic sources and the measured responses. The corresponding mathematical problem is ill-posed as is usual for inverse problems [40, 107]. The goal of this thesis is to derive fast solvers for time-harmonic Maxwell's equations for low or moderate frequencies. Toward this goal, key hindrances in the forward problem are identified that point to analytic reformulations of the partial differential equations (PDEs) in terms of scalar and vector potentials. This analytic model lends itself to a straightforward discretization applying finite-volume techniques. The discretization has conservative properties on the discrete level typical of finite-volume methods and also successfully approxi-mates interface conditions imposed by discontinuities in the material param-eters. Finally, a set of preconditioners for Krylov-subspace methods applied to the resulting non-Hermitian sparse linear systems of algebraic equations is explored. Of particular interest is a multigrid preconditioner that exploits both the structure of the P D E problem at hand and well-established solvers for elliptic P D E problems [19, 35, 61, 60, 92]. The combination of the afore-mentioned efforts leads to a robust solver for the forward-modeling equations. Although the analysis and techniques presented here focus explicitly on geo-physical applications, extensions exist to other low-frequency electromagnetic applications such as medical imaging. 2 1.1 Background Sketch The past decades have seen the development of many numerical tech-niques for the computation of electromagnetic fields as well as significant ad-vances in computer hardware, both in processor speed and memory capacity. Fast and accurate methods for predicting electromagnetic phenomena are more readily available. Three-dimensional numerical simulation of the interaction of electromagnetic fields with matter is attainable on modern desktop machines. The basic physics of electromagnetic fields is governed by Maxwell's equations. Among alternative mathematical formulations of Maxwell's equa-tions, integral equations (IEs) and partial differential equations (PDEs) lead to very different computational approaches. With computational methods derived from IEs, a three-dimensional boundary-value problem reduces to a two-dimensional problem over the boundary of the domain of interest (e.g. boundary element methods or the method of moments [59, 98]). However, even with a significant reduction in the number of unknowns, the compu-tational cost of generating the full system matrix and difficulties in solving the linear equations often makes this approach more costly than comparable P D E methods [77]. It is also difficult to formulate the appropriate IE for geo-metrically complex inhomogeneous structures, possibly requiring a nontrivial derivation of a geometry-dependent Green's function; P D E solutions are much more flexible in dealing with such media. Thus, P D E formulations of electro-magnetics and consequent numerical techniques are the focus of this thesis. For static problems, scalar and vector potentials for Maxwell's equations 3 are commonly used [69]. As observed in [16], the paradigmatic equation for electromagnetics in the static case was long thought to be the diffusion problem div (e grad <p) = p for the electrostatic potential 0. During the 1970's, advances in finite ele-ment methods popularized simulation and analysis for diffusion problems of this form [23, 69, 102]. In two dimensions, the corresponding equation for magnetostatics curl (//_1curl a) = J has only one component for the vector potential a and reduces to a similar scalar div-grad problem. Thus, it was commonly perceived that the div-grad and the curl-curl equations are effectively equivalent. However, in three dimensions, this is not the case; complications including the nontrivial kernel of the curl operator and coupling of field components in three dimen-sions make the latter operator much more difficult to handle. This realization led to intense research activity into suitable computational methods in three dimensions (e.g. edge elements [11, 16]). Although finite-element methods are rich in theoretical tools and pro-vide great geometric flexibility, finite-difference methods and finite-volume methods are simpler to implement and can still provide accurate descriptions of solutions for many realistic electromagnetic phenomena [6, 25, 44, 79]. The first finite-difference codes to model Maxwell's equations in three dimensions were based on the well-known Yee method [115]. This spatial discretization 4 scheme uses a staggered grid to explicitly enforce discrete analogues of con-servation properties; this is important for certain applications including elec-tromagnetic scattering, radiation, and waveguides [69, 104]. Although used mostly in electrical engineering applications, Yee's discretization has over the last ten years been applied in both time- and frequency-domain problems by geophysical researchers [3, 84, 100, 109]. Beyond the need for accurate discretizations in electromagnetic simula-tion, the large sparse systems of linear equations that typically arise from the discretization of PDEs in three dimensions must be solved efficiently. Direct methods have long been known to be inadequate with the exception of certain cases with unique structure (see, e.g., [48]). With the dramatic increases in processor speeds and growth of memory capacity in modern computers, com-putations on a large scale are more common and iterative methods become crucial [10, 49, 93]. Although originally conceived as a direct method, the widely-known method of conjugate gradients has had a great impact as an iterative method for the solution of linear systems obtained from discretizations of elliptic P D E problems [60, 61]. One appealing feature is that the convergence behavior of the conjugate-gradient method for elliptic problems is understood. Specifi-cally, for a square linear system Hx = b with a Hermitian positive-definite matrix H, the number of iterations required for the residual b — Hxi of the approximate solution 2 ; to converge within a set tolerance is 0(y/cond(H)) = 0 ( v / A m a x / A m i n ) , where cond(jrY) is the ^-condition number of H and A m a x and 5 A m i n are respectively the largest and smallest eigenvalues of H [49, 93]. For discretizations of second-order elliptic PDEs with grid spacing h, the bound on the number of iterations typically increases as 0(h~l) since the condition number is 0(h~2) (see, e.g., [49, Chap. 9]). Unfortunately, a great number of linear problems are non-Hermitian and indefinite or both, including discrete systems associated with Maxwell's equations. Throughout the 1980's and 1990's, much effort went into the de-velopment of Krylov-sub space methods for the solution of more general lin-ear algebraic equations. Methods such as GMRES, QMR, and BiCGStab [10, 45, 49, 93, 94] provide generalizations of the conjugate-gradient method that apply for non-Hermitian problems. However, the theory guaranteeing convergence behavior of non-Hermitian matrix iterations is scarce. In spite of the lack of such theoretical results, non-Hermitian Krylov-subspace itera-tions are applied in most areas of physical science with great success (see, e.g., [6, 26, 28, 54] for examples relating to electromagnetics). In many ways, more important than the search for more variants of Krylov-subspace methods themselves is the search for suitable preconditioners for matrix iterations. Given an invertible square matrix A, a preconditioner is a matrix M that in some way approximates A"1 so that MA (or AM) has better spectral properties than A itself. The preconditioner M is selected so that the matrix-vector product Mz is cheap to evaluate for any column vector z (e.g. by F F T , fast Poisson solvers, ADI methods, multigrid methods, etc.) since that product has to be computed within each iteration of a Krylov-subspace 6 method. While preconditioners based on classical stationary iterative methods are simple to implement in a black-box solver and may have other advantages, many important preconditioners exploit properties of the continuous problem underlying the discrete linear system. For comprehensive summaries of Krylov-subspace methods and some standard preconditioners, see [10, 49, 93]. Over the last thirty years, an important class of iterative methods for elliptic P D E problems emerged, namely the family of multigrid methods [17, 18, 19, 22, 56, 92]. Multigrid methods rely on the fact that the error of a fine-grid discrete approximation of the solution of an elliptic P D E can be split into high- and low-wavenumber components; the high-wavenumber com-ponents are damped rapidly by the action of a smoother (e.g. iterations of a classical stationary iterative method like Gauss-Seidel or SOR) on a fine grid, while the low-wavenumber components can be eliminated on a coarser grid where the number of discrete unknowns is fewer. Given a discretization of an elliptic PDE, the key components of a multigrid method are a hierarchy of fine and coarse grids, a basic smoothing or relaxation method for improving discrete solution approximations, a re-striction operator and an interpolation operator for transferring grid functions between grids, and a coarse-grid operator approximating the action of the dis-crete operator on the coarse grid (see [18, 22] for more details). A multigrid method is constructed based on the following two-grid heuristic: 1. Pre-smoothing: Given a guess for the discrete solution resolved on the fine grid, smooth out the error with a few sweeps of the smoother. 7 2. Restriction: Compute the residual of the current solution on the fine grid and project the result onto the coarse grid with the restriction operator. 3. Coarse-grid solve: Obtain the coarse-grid error by applying the inverse of the coarse-grid operator to the coarse-grid residual. 4. Prolongation: Interpolate the coarse-grid error up to the fine grid with the interpolation operator, giving an approximation of the fine-grid error of the current fine-grid solution. 5. Post-smoothing: Update the current fine-grid solution using the fine-grid error and perform a few more smoothing sweeps to improve the solution. At the level of the coarsest grid, the appropriate grid equations are solved exactly or nearly exactly to eliminate the low-wavenumber components of the error. A multigrid method is based on a recursive implementation of the above two-grid scheme on a hierarchy of grids (i.e. applying the same two-grid heuris-tic for the coarse-grid solve in step 3). Such a recursive iteration is called a V-cycle since problems are solved on progressively coarser grids until reaching the coarsest grid, after which the grids are revisited in reverse order. A vari-ant of a V-cycle where the coarsest grid is visited twice in a single iteration is called a W-cycle. Multigrid methods stand out amongst other iterative methods applied to discretizations of elliptic PDEs since well-constructed multigrid methods have grid-independent rates of convergence. Specifically, for a problem dis-8 cretized on a fine grid of spacing h, if a multigrid method is constructed and applied, then the number of multigrid W-cycles required to reduce the error within a set tolerance is constant as the grid spacing h is decreased [17, 19]. The same is not true for Krylov-subspace methods, such as the method of conjugate gradients, where the corresponding number of iterations increases as 0(h~l) for standard discretizations of elliptic PDEs [10, 49]. Thus, although multigrid software can be tricky to implement, for sufficiently large problems (i.e. problems resolved with small grid spacings), the pay-off can be highly significant. For more detailed surveys and background on multigrid methods, consult [17, 18, 19, 22, 56, 92]. 1.2 Outline of research A background survey of electromagnetic theory needed to describe ba-sic electromagnetic phenomena is presented in Chapter 2 with a particular focus on geophysical problems. Starting from Maxwell's equations, various physical assumptions used in practice are described, including, among others, constitutive relations, time-harmonic and quasistatic models, and boundary conditions. A clear understanding of the modeling assumptions imposed is crucial from physical, mathematical, and computational perspectives; it is de-sirable that all assumptions yield a description that is consistent with the ac-tual physics, that is mathematically sound and solvable, and that is amenable to reliable computational techniques. Thus, in Chapter 3, the forward-modeling problem is formulated in the 9 frequency-domain drawing particular attention to the physical properties of the systems being modeled. First, the frequencies of interest are of moderate magnitude, so the effective wavelengths of the scattered fields are much greater than the length scales considered. Secondly, the material properties to be recovered in the inverse problem are isotropic and inhomogeneous with jump discontinuities of several orders of magnitude occurring at media interfaces. Finally, the conductivity coefficient effectively vanishes in the air, so traditional P D E formulations, although well-posed under suitable assumptions, lead to ill-conditioned linear systems when discretized with standard techniques. Although the latter two of the preceding observations are hindrances, the first is actually an advantage. For sufficiently low frequencies (and param-eter ranges suitable for geophysical problems), the prevalent flow of solutions of Maxwell's equations in the time-domain is parabolic rather than hyperbolic; that is, the physics of the electromagnetic field is prevalently diffusive rather than dispersive as in high-frequency electromagnetics. This feature is eventu-ally exploited to construct preconditioners, but it relates also to the selection of suitable analytic formulations of Maxwell's equations. Scalar and vector potentials are introduced with a Coulomb gauge con-dition in Chapter 3 to circumvent the problems induced by the dominant curl-curl operator in the frequency-domain P D E . The resulting indefinite system is then stabilized by addition of a vanishing term that lies in the kernel of the dominant curl operator. Finally, an extra differentiation recasts the P D E system in a diagonally-dominant form reminiscent of a pressure-Poisson 10 formulation for incompressible fluid flow. Unlike the related problem in fluid flow, there is some freedom to choose boundary conditions for the potential fields so that this formulation leads to robust discretizations. Sections 3.2 and 3.3 fill in the remaining details by deriving suitable boundary and inter-face conditions for the scalar and vector potential fields and by providing a discussion of a related weak formulation. After recalling the Yee (or FD-FD) discretization at the beginning of Chapter 4, a finite-volume method is derived with discrete scalar and vector potential fields as unknowns on a staggered grid [53, 54]. An important feature of the discretization is that it retains the important conservation properties of Yee's discretization while pointing to logical ways to represent the material inhomogeneities within the discrete equations. In particular, the conductivity at the interface of discontinuity between adjacent cells is calculated using har-monic averages, while the permeability at the edges between adjacent cells is given by arithmetic averages. The discretization of material parameters typically used throughout the geophysical community consists of a grid composed of cells within each of which the material parameters are assumed constant (e.g., [37, 84]). While convenient for inverse problems, this model can only resolve interfaces be-tween distinct media with a staircase approximation of the actual geometry. For high-frequency problems, failure to resolve the geometry accurately may cause great difficulty in practice (e.g. [105]). However, for low-frequency prob-lems, the diffusive nature of the P D E operator tends to smooth out computed 11 errors and the solutions are not as sensitive to perturbations in the coefficients. In Section 4.4, some analysis of a simple inhomogeneous diffusion problem is presented to estimate the nature of the bounds on the perturbations in the fields induced by staircasing. The analysis is accompanied by supporting ex-periments in three dimensions. In Chapter 5, a new family of preconditioners based on sparse approx-imations of the original linear systems is constructed. Specifically, in the low-frequency limit, the diagonal blocks of the linear system are dominant and thus suggest an obvious block diagonal preconditioner. The diagonal blocks are effectively discretizations of scalar Poisson and inhomogeneous diffusion problems in three dimensions. The blocks can be inverted exactly, or approx-imate inversion can be achieved either by Incomplete-LU (ILU) factorization or by a multigrid solver such as B O X M G [35, 36]. All of these strategies are compared in the experiments presented. Since the coefficient matrix is complex and non-Hermitian, the perfor-mance of Krylov-subspace methods on the preconditioned system is not well understood. However, a Fourier analysis yields frequency-dependent bounds on the spectrum of the preconditioned system from which a condition bound can also be derived under some simplifying assumptions [5]. The analysis, although limited by the assumptions, does indeed capture the local behavior of the discrete system being solved. Indeed, the numerical tests show grid-independent rates of convergence of the resulting multigrid-preconditioned solver in spite of the non-Hermitian character of the linear system. The ef-12 fectiveness of analyzing local behavior of the discrete operators in this case follows largely due the physical assumptions inherent in the model. Finally, a summary of this research and some future research directions are presented in Chapter 6. 13 Chapter 2 Background Theory of Electromagnet ism James Clerk Maxwell elegantly formulated the classical theory of elec-tromagnetism in 1864. Maxwell unified previous theoretical and experimental knowledge of electromagnetic phenomena into a more general and advanced theory summarized in a system of partial differential equations (PDEs) known as Maxwell's equations. Physical scientists have studied this particular model of macroscopic electromagnetic phenomena extensively for over a hundred years, laying out the theoretical foundation on which modern electronic tech-nology is built. 2.1 Maxwell's Equations Maxwell's equations can be formulated using the language of differen-tial forms [16, 67, 73], but the more common language of vector analysis suf-fices. Written as partial differential equations in Minkowski form (see [68, 76]), 14 Maxwell's equations are curl E + dtB = 0, (Faraday's law) (2.1a) curl i f - dtD = J, (Maxwell-Ampere law) (2.1b) div D = pe, (2.1c) d i v B = 0. (2.1d) The quantities represented in (2.1) are distinct fields that are functions of spa-tial position x and time t. The vector fields are E (the electric field intensity), D (the electric flux density), H (the magnetic field intensity), B (the mag-netic flux density), J (the electric current density), and pe (the electric charge density). The four fields E,D,H, and B describe the total electromagnetic field. The PDEs (2.1) hold in an inertial frame in which the material media are at rest and the charges that interact with the electromagnetic field can move and produce currents [68, 101]. Table 2.1 summarizes these fields in terms of fundamental dimensions of mass (M), length (L), time (T), and current (I). Quantity Units ~ E Electric field intensity volts/meter (V/m) M L T " 3 / - 1 D Electric flux density coulombs/meter2 (C/m 2 ) L~2TI H Magnetic field intensity amperes/meter (A/m) L - 1 / B Magnetic flux density webers/meter2 (W/m 2 ) MT~2I~l J Electric current density amperes/meter2 (A/m 2 ) L~2I pe Electric charge density coulombs/meter3 (C/m 3 ) L~3TI Table 2.1: Field quantities in Maxwell's equations. There is an additional continuity equation d i v J + d tp e = 0 (2.2) 15 that expresses the conservation of electric charge. In three dimensions, only 7 of the 9 scalar equations (2.1) and (2.2) are independent since (2.1d) follows directly from (2.1a) (with the assumption divJB = 0 at t = 0) and (2.2) is a consequence of (2.1b,c). Usually, (2.1a-c) or (2.1a,b) and (2.2) are chosen as the independent equations. 2.2 Derived Models Maxwell's equations (2.1) in their most general form describe a great range of phenomena. There are a variety of physical models in applied sciences that arise from some simplifying assumptions in (2.1). Among possible models, there are a number of further assumptions that can be made: 1. Constitutive relations 2. Stationary models: electrostatics and magnetostatics 3. The time-harmonic model 4. The quasistatic model 5. Electric source currents 6. Magnetic source currents 7. The magnetotelluric model Of these possible modeling assumptions, the constitutive relations (2.3) are crucial. The time-harmonic model is considered in this thesis, occasionally 16 making use of the quasistatic assumption. Stationary models are not studied in this thesis, although much of the standard theory developed for static problems (in particular, the use of potentials) is applied. Assumptions dealing with boundary conditions are discussed in Section 2.3. 2.2.1 Constitutive Relations Since there are 7 independent scalar equations from (2.1-2.2) that involve 16 scalar unknowns (including the components of the vector fields), the system is under deter mined. A determinate system requires further assumptions. To-wards this end, impose constitutive relations between the field quantities in order to make the system (2.1) definite [16, 76]. These take the form D = eE, (2.3a) B = pH, and (2.3b) J = oE + J e , (2.3c) where e is the electric permittivity, p is the magnetic permeability, and cr is the specific conductivity of the macroscopic media in which the electromagnetic field exists. The current density J is decomposed into «7e (the current density due to some external applied electric source), and oE (the current density induced in conducting matter by the source current J e ) . The quantities in (2.3) are summarized in fundamental units of mass, length, time, and current in Table 2.2. The quantities e, and o are macroscopic properties of matter. They are generally tensor fields that depend on the fields E and H as well as the 17 Quantity Units ~a Specific conductivity siemens/meter M _ 1 L _ 3 T 3 / ' 2 e Electric permittivity farads/meter M~lL~zTAI2 H Magnetic permeability henrys/meter MLT~2I~2 Table 2.2: Constitutive parameters for Maxwell's equations. spatial variables x and the time t. The tensor fields e and fj, are positive def-inite while the tensor field a is positive semi-definite (because it is physically plausible for a to vanish). For isotropic media, the tensor fields for e, fi, and o reduce to scalar multiples of an identity tensor, so they are essentially scalar fields; for linear media, they are independent of the state of the electromag-netic field (i.e. of E and H); for homogeneous media, they are independent of position x; and, for dispersive media, they depend on the frequencies of the electromagnetic fields. In free space, e and // are isotropic and homogeneous; the corresponding permittivity of free space is denoted eo and has the value e0 := 8.85 x 10~12 F / m in SI units, while the permeability of free space is denoted /J,O and po := 47T x 10 - 7 H/m. Free space is nonconductive, so o = 0 in free space; El-ementary analysis of propagation of electromagnetic waves within isotropic, homogeneous media allows exact solutions of Maxwell's equations to be de-termined [101]. For more realistic physical simulations of electronic devices, for scattering experiments, or for geophysical data inversion, the material ten-sors are at the very least inhomogeneous, usually with discontinuities across material boundaries. It is typical to choose either (E, H) or (D, B) as the unknown fields 18 once the constitutive relations (2.3) are assumed. Opting for the unknown fields (E,H), Maxwell's equations (2.1) become curlE + pdtH = 0, (2.4a) curl H - o E - edtE = Je, (2.4b) div (eE) = pe, (2.4c) div (pH) = 0. (2.4d) The system (2.4) is a linear, first-order, hyperbolic system of PDEs [89]. Within this work, the material properties e, p, and a are assumed to be in-dependent of time t and of the frequencies of the electromagnetic fields (i.e. nondispersive); this permits moving the quantities e and p through the oper-ator dt in (2.4a,b) without concern. However, great care is required in manip-ulating the operator div in (2.4c,d) unless e and p are spatially constant (i.e. homogeneous). 2.2.2 Stationary Models: Electrostatics and Magneto-statics The most obvious simplification is to assume that the electromagnetic fields do not vary in time. As such, all terms with time derivatives vanish in (2.4). An immediate consequence is the electrostatic model curl E = 0, D = eE, (2-5) divZ? = pe. 19 Assuming the static charge density is known, it is possible to identify the elec-tric field E as the gradient of a scalar potential field, so E = — grad</>. With this notation, the essential model of linear electrostatics is the inhomogeneous diffusion equation -d iv (e grad cp) = pe. (2.6) With a suitable domain and boundary conditions specified, the diffusion P D E (2.6) yields a well-posed boundary-value problem [46, 89]. In the case where e = eo (i.e. the permittivity is everywhere constant) and where the domain is all of M3, standard analysis (e.g. [101]) gives the electrostatic potential 47re 0 JJJus \x - y\ In a similar way, the model of magneto statics is given by curl H = J, B = fiH, (2-8) divB = 0, for some known time-independent source current density J (the effects of conduction are not included in this model, so there is no oE term). The vanishing divergence of B usually implies the existence of a vector potential field for B, so H — / / _ 1 curla for some vector field a. Thus, the essential equation describing linear magnetostatics is the double-curl equation curl (Ai _ 1 c u r l a ) = J . (2.9) 20 If p = po, the vector potential over all space is The magnetic field can be recovered "naturally with the Biot-Savart law [43, 76] 2.2.3 The Time-Harmonic Model For this model, all the scalar fields and the components of the vector fields in (2.4) are assumed to be of the form u(x,t) = Re(u(x,uj)e~lU}t) for some constant frequency UJ G R. The resulting model of electromagnetism is commonly called time-harmonic. Then, Maxwell's equations (2.4) become curl E - lupH = 0, (2.12a) c u r l H - a \ E = J e , (2.12b) div (eE) = pe, (2.12c) div (pH) = 0. (2.12d) In (2.12b), the complex conductivity o := a — luie (2-13) is introduced for notational convenience. The equations (2.4) are Maxwell's equations in the time domain, whereas the equations (2.12) are Maxwell's equations in the frequency domain (sometimes called time-harmonic Maxwell's 21 equations). One advantage of the frequency-domain formulation of Maxwell's equations is that, in principle, solutions of (2.12) can be found for a few key frequencies of interest. Formally, the system (2.12) is obtained by substituting —itue for every occurrence of dt in (2.4). Within this thesis, the same notation is used for fields u = u(x,t) in the time domain and the corresponding fields u = u(x, uj) in the frequency domain; the meaning is always clear from the context. 2.2.4 The Quasistatic Mode l For many physical applications, it is possible to neglect the Maxwellian displacement current term IUJD from (2.12b) (or dTD from (2.1b) in the time-domain formulation); in the frequency-domain, this amounts to substituting cr <— a in (2.12b). This quasistatic assumption (also called the eddy-current ap-proximation) is reasonable assuming that the displacement current IUJD is neg-ligible relative to the other terms in the Maxwell-Ampere law (2.12b). One key indicator of the validity of the quasistatic assumption in the frequency-domain model is the magnitude of the ratio u>e/o; for low-frequency experiments within conductors, this ratio can vary from 10~16 to 10 _ 1 (see [16, 111]). In many high-frequency electromagnetic applications (e.g. waveguides) or applications in vacuum (e.g. radar scattering), the displacement current's contribution to the evolution of the electromagnetic field cannot be ignored as readily. For a discussion of parameter ranges in which the quasistatic assumption is valid, see [16, 111]. 22 The mathematical importance of the quasistatic assumption is that it changes the underlying character of Maxwell's equations. Ordinarily, the time-domain equations (2.4) constitute a hyperbolic system of PDEs; dropping the displacement current dtD makes the P D E system parabolic. The ratio of the magnitude of the conduction current aE to the magnitude of the displacement current dtD — edtE reflects the relative importance of diffusion to dispersion in transmission of electromagnetic energy on the time scales considered. The balance between these physical processes strongly influences the development of mathematical models and discretization schemes in low- and high-frequency electromagnetics. 2.2.5 Electric Source Currents Generally, the conduction currents within matter as described by Max-well's equations (2.4) (2.12) are induced by some time-varying field driven by some external electric current density Je. Within the domain of geophysical applications, certain source currents are naturally occurring, such as iono-spheric currents that cause geomagnetic variation, while others are artificial, such as current loops used in geophysical prospecting [111]. In either case, the source of electromagnetic induction is incorporated into Maxwell's equations through the constitutive relation (2.3c), namely J = aE + Je. In the following, let ex, ey, and ez denote the orthogonal Cartesian coordinate unit vectors. Let x := (x, y, z)T denote the coordinates of a point 23 in space. The function 5 refers to the usual Dirac-<5 function and H refers to the Heaviside-step function (see e.g. [101]). The sources are understood as distributions that act on a space of testing functions (see [46, 89] for definitions of distributions). Examples of source fields include: 1. A vertical electric dipole Je(x,t) = I(t)5(x)S(y)(H(z -£)- H(z))ez. (2.14a) In (2.14a), £ > 0 is the separation of the electrodes between which a potential difference is maintained that generates the dipole field. The function I(t) describes the current that flows between the electrodes as a function of time. 2. A vertical magnetic dipole Je(x,t) = I(t)8(z)[(H(x)-H(x-£))(5(y)-5(y-£))ex - (8(x) - 5(x - £)) (H(y) - H(y - £)) ey]. (2.14b) The magnetic dipole field is generated by a current loop with normal direction ez. Again, I(t) is the current flowing through the loop as a function of time and £ is the sidelength of the loop. 3. A source current along a line Je(x,t) = I(t)5(y)5(z)ex. (2.14c) 24 The current I(t) flows in the x-direction. Notice that, unlike the previous source fields, this field has unbounded support. Other source fields of unbounded support include current sheets in a horizontal or vertical plane [111]. The sources in (2.14) are used in time-domain modeling. If the current inducing the source is of the form I(t) = i oe _ t / V in (2.14b-d), the sources are referred to as transient (that are of interest in geophysical prospecting [110]). For frequency-domain models, the applied currents are of the form I(t) = Ioe~llJJt in (2.14b-d). Exact solutions of Maxwell's equations are known for all of the above with time-harmonic sources in vacuum and in half-space conductivity models due to the high degree of symmetry [76, 101, 110, 111]. 2.2.6 Magnetic Source Currents Given the apparent symmetries in the curl equations (2.12a-b), it is reasonable from a mathematical perspective to place some kind of a source term on the right-hand side of (2.12a). Under this assumption, the time-domain curl equations become curl E + pdtH = J, (2.15a) curl If - (a + edt)E = J, (2.15b) div (eE) = p( (2.15c) div (pH) = p, (2.15d) 25 and in the frequency domain, curl E — tupH = J , (2.16a) curl H - aE = J. (2.16b) div (eE) = p, (2.16c) div (pH) = p, (2.16d) In (2.15) and (2.16), Jm is a magnetic current density and pm is a magnetic charge density. Such magnetic currents and charges are not physically realistic; magnetic charges correspond to magnetic monopoles that have never been observed experimentally. As a result, J m and pm are null in all physically relevant situations [43, 68, 101]. However, (2.15) and (2.16) are useful for the purpose of analysis and derivation of numerical schemes. For instance, in the forward-modeling prob-lem, if there is a known primary field (Ep, Hp) in the background, it may be more convenient to solve for the unknown secondary field (Es, Hs), assuming that the total electromagnetic field is given by (Ep + Es, Hp + Hs) [6, 14]. Solving for a secondary field may be more useful in the event that the secondary field is significantly smaller in magnitude than the primary field. Magnetotel-luric models provide an example of a physical model in which secondary fields and hence magnetic source currents enter the mathematical formulation of Maxwell's equations. 2.2.7 The Magnetotelluric Model 26 In a large class of geophysical problems, the flow of electromagnetic fields in the higher atmosphere induces currents in conductive matter under-ground; this is the basic principle of the magnetotelluric (or MT) experiment (see [71, 77, 95]). To derive the basic physical model for the M T experiment, assume that the total electromagnetic field satisfies Maxwell's curl equations in the frequency domain with no applied source currents, i.e., c u r l E - zujfiH = 0, (2.17a) c u r l H - ( a - zue)E = 0. (2.17b) The equations (2.17) are assumed to hold throughout all of space. The material tensors cr, e, and p are isotropic, but inhomogeneous; in particular, they are assumed to be piecewise constant scalar fields throughout M 3 . Assume that the field admits the decomposition E = Ep + E\ and (2.18a) H = Hp + Hs. (2.18b) The field (EP,HP) in (2.18) is a primary field that is assumed to satisfy Maxwell's equations with some known coefficients and source terms. If the primary field is known, then solving the original system (2.17) amounts to solving a similar P D E for the secondary field (Es, Hs) induced by the primary field. The left-hand side of the P D E for the secondary field is determined by the primary field and describes the response of the secondary field due to the primary one [6, 14]. For the magnetotelluric modeling, it is sensible to use a half-space or a • 27 1-D layered-earth solution (i.e., piecewise constant functions of z only, see [111, 110]) of Maxwell's equations for the primary field based on the assumption that the earth being modeled has that basic structure far away from the region of interest [71]. Thus, assuming that the primary field satisfies Maxwell's equations with no source currents in a half-space model, then curl EP = 2UJ/1HP, (2.19a) cu r l Hp = (a- ILO7)EP (2.19b) where a, e, and JI are models for a half-space or ID layered earth. The equa-tions (2.19) admit analytic solutions provided that the coefficients in (2.19) have such a simple one-dimensional structure. Under the assumption (2.19), since J e = 0 for the primary field, the secondary field satisfies cur l ES - iufiHs = J m , (2.20a) cur l H" - (a - icue)ES = J e , (2.20b) where J m := IU (p - ji) Hp, (2.21a) J e := [(a - a ) - i u ( e - e)] E P . (2.21b) In decomposing the total field, the P D E for the secondary field has a mag-netic current source term J m as well. The sources in (2.21) generally have unbounded support since the fields EP and Hp will generally have unbounded support. However, if it is known a priori that o = o, e = ?, and fx = Ji outside 28 some bounded domain, then the source terms indeed have compact support. In fact, if the computational domain is large enough, the normal and tangen-tial components of J m and Je in (2.21) on the boundary of the computational domain fl vanish, and these fields do not affect the boundary conditions. If there is a more complicated structure to the material tensors outside fl, then the source terms may enter into the boundary conditions. Suitable boundary conditions are presented in the next section. 2.3 Boundary Conditions The forward-modeling problem refers to determining the unknown fields E, H, and pe under the assumption that the quantities u> (in the frequency domain), J e , e, p, and a are all known. The forward-modeling problem can be solved in the time domain with (2.4) or in the frequency domain with (2.12). In either case, the system consists of 8 differential equations involving 7 un-known quantities. In practice, the divergence equation (2.4c) (or (2.12c) in the frequency domain) defines the charge density pe once the field E is deter-mined. Moreover, (2.12d) can be derived from (2.12a) under the assumption UJ 7^ 0. Thus, consider (2. ld) as a differential invariant and use the curl equa-tions (2.4a,b) (or (2.12a,b)) as the system of PDEs determining E and H in the forward-modeling problem. Even with the constitutive relations (2.3) and various modeling assump-tions about the source fields, a suitable set of boundary and initial (in the time domain) conditions is needed to make the forward-modeling problem 29 well posed. For an unbounded domain ficR3, suitable boundary conditions at infinity are well known [69, 91]. However, the question of what boundary conditions to apply on the boundary dfl of a bounded computational domain O C I 3 is actively debated due to the inherent difficulty of balancing the need for physically relevant boundary data with the desire for an elegant, consistent, mathematical formulation. This direction of research in electromagnetics par-allels the search for open boundary conditions in computational fluid mechan-ics [51, 50]. Absorbing boundary conditions have been used for over twenty years in electromagnetic codes to absorb outgoing waves at dfl and to ensure that reflected waves do not corrupt the solution [41]. As a first-order approxi-mation, absorbing boundary conditions correspond to Sommerfeld's radiation condition for the Helmholtz equation; at high frequencies where greater ac-curacy may be required, higher-order absorbing boundary condition are also possible [58, 57, 69]. A different class of approximate boundary conditions stems from a perfectly matched layer technique (PML) [12] in which some absorbing matter is placed immediately outside the artificial boundary dfl to absorb outgoing waves; for a numerical study of implementations of various artificial boundary conditions, see [104]. It is not necessary to assign boundary conditions on the artificial bound-ary with certain formulations. One alternative is to use boundary integrals (e.g. boundary element methods [59]) to connect the problem within the finite computational domain to the unbounded exterior region. The use of infinite elements also allows a computational framework that incorporates the behav-30 ior of the fields at infinity into the solution [47]. Both of these techniques have some advantages but may introduce other computational difficulties. For a detailed review of different techniques used for the numerical solution of problems on unbounded domains, see [106]. It is important to recognize that for the classes of geophysical applica-tions studied herein, the time-harmonic fields are driven by low- or moderate-frequency sources; as such, the influence of boundary conditions at the ar-tificial boundary dQ, is minimal. This is in contrast with the situation in high-frequency scattering applications where errors in the computation of the electromagnetic field at dfl contaminate the entire solution. For the diffusion-dominated flow of the electromagnetic field on geophysical length scales, the accuracy of the artificial boundary condition is not as crucial as it would be at higher frequencies. This thesis focuses on frequency-domain formulations of Maxwell's equa-tions. As such, the discussion that follows is limited to well-posed P D E prob-lems in the frequency domain. For results describing initial and boundary conditions for well-posed P D E problems based on Maxwell's equations in the time domain, see [4, 11, 70, 73, 97]. 2.3.1 Frequency-Domain Boundary Conditions First, consider the time-harmonic P D E problem (2.12) specified in an 31 unbounded spatial domain flcK3. Then, as |cc| —> oo, provided that \x\E(x,ui) and \X\H(X,LO) are bounded and that ' ' V ' 1 ' V (2.22) i r i / \ , XXH(X,OJ) / i \ E(x,u) + \a^2 = o{W\) uniformly in all directions, then the frequency-domain problem (2.12) is guar-anteed to have a unique solution [97]. The conditions (2.22) are the Silver-Miiller radiation conditions that are analogous to the Sommerfeld radiation condition for the Helmholtz equation [69, 91]. On a bounded domain f l c R 3 , suppose that the boundary dfl is par-titioned into two disjoint parts, dflE and dVtH := dfl\dflE, one of which may be empty. Then, the following inhomogeneous boundary-value problem is well-posed (see [16, chap. 9] or [15]): curl E — lupH — Jm ( i e (1); curl H - oE = Je (xeSl); V 1 (2.23) n x E = n x e (x e d£lE); n x H = n x h (x E dflH), where e and h are some known fields whose tangential traces are well-defined. The explicit spaces in which the solution (E, H) is sought are detailed in Section 2.5 where equivalent integral formulations are specified. The specification of the tangential trace n x E = 0 is the boundary condition that applies at the boundary of a perfect electric conductor (PEC); hence, this boundary condition is often called a P E C boundary condition. The corresponding boundary condition n x H = 0 is the boundary condition for a perfect magnetic conductor (PMC). While both P E C and P M C boundary 32 conditions represent idealizations that do not exist in nature, they are reason-able approximations and are used in many electromagnetic models for their simplicity [11, 63]. A standard first-order absorbing boundary condition is of the form where g is the tangential trace of some known field [41, 97]. The particular absorbing boundary condition (2.24) ensures that electromagnetic waves nor-mally incident on the boundary will be completely absorbed. This is more physically reasonable than a P E C or P M C , although there will be some reflec-tions due to waves striking the boundary obliquely. However, if the wavevector k of an incident electromagnetic wave is known (for instance, the wave that generates the primary solution in the M T experiment in Section 2.2.7), a slightly different absorbing boundary condition applies, namely - n x n x ( e - — ) E + A(k,n)n x yAH = g (xedfl), (2.25) where A(k, n) is a symmetric positive-definite matrix chosen so that waves on 80, with wavevector k are completely absorbed [97]. Well-posedness of time-harmonic Maxwell's equations with boundary conditions (2.24) or (2.25) is shown in [95, 97]. Higher-order absorbing boundary conditions are discussed in [69]. 2.3.2 Interface Conditions n x (—n x E + H) = n x g {x <= dfl), (2.24) 33 The material tensors e, / i , and a are generally discontinuous at interfaces between distinct kinds of matter. Due to such discontinuities, the differential equations (2.4) and (2.12) (and their generalizations (2.15) and (2.16)) do not apply at the interfaces; they hold strictly within regions wherein the material tensors are continuous. Thus, in addition to boundary conditions that de-scribe the behavior of the electromagnetic field at the artificial boundary 30, interface conditions apply at material interfaces. To derive the associated interface conditions, let T C fl be any surface that divides fl into two disconnected portions, O 1 and O n . Consider an in-finitesimal cylinder Vs C fl of height 5 > 0 that straddles T. The cylinder Vs intersects T in a surface IV Let the boundary of Vs consist of three disjoint parts: the sides, top, and bottom respectively (see Figure 2.1), so 8V5 = dV6sides U dV} U dVsn. (2.26) Finally, let n 1 and n 1 1 denote the unit normal vectors on dV5l and dVj1 respec-tively with the conventional outward orientation. Vs Figure 2.1: Infinitesimal cylinder Vs for deriving interface conditions. Multiplying the divergence equation (2.1c) by some arbitrary smooth 34 function ip and integrating by parts yields / ibnl-DdS+ f ipnu-DdS= [ tppedV Jav} Jdv}1 Jv5 + [ gradip-DdV- f ipn-DdS. (2.27a) Jvs JdVf'dcs Choose the test function ip in (2.27a) so that tb = 0 on dVssides. Then, / tbnl-DdS+ f ipnu-DdS= [ ippedV JdV} JdVp Jvs + f gradib-DdV. (2.27b) Jvs As 5 I 0, the second volume integral on the right-hand side of (2.27b) vanishes assuming that D and g r a d ^ do not have surface singularities. The other volume integral on the right-hand side reduces to a surface integral of a surface electric charge density. The preceding argument holds for any surface Y C 0 (not just interfaces), so it follows that n 1 • [Dl - Du] = ge, (2.28a) holds pointwise on T, where Qe is the surface electric charge density. The expression Dl in (2.28a) is the limiting value of D approaching a point on T from within Q1 and Du is similarly defined. In a similar manner (see, e.g. [68, 69, 76, 101]), the conditions n 1 • [B l - Bn] = Qn, (2.28b) n I . [J 1 - J«] = -dtQe, (2-28C) n 1 x [Hl - H11} = j e > (2.28d) n 1 x [El - E11] = Jm (2-28e) 35 also apply across any surface T C 0. In (2.28), ge and gm are surface elec-tric and magnetic charge densities while je and jm are surface electric and magnetic current densities respectively. The interface conditions (2.28) are derived explicitly assuming the fields E, H, D, and B do not have singularities in any of their components. This modeling assumption rules out the possibility of singular layers that might occur in cracks or tiny air gaps in a conductor (see [16]). For the geophysical problems considered in this thesis, this assumption is reasonable. Further, even where magnetic charges and currents are formally used for mathematical purposes (cf. Section 2.2.6), the surface densities gm and jm are assumed to vanish, assuring the continuity of n • B and n x E throughout fl. Moreover, between media of finite conductivities, no surface electric current density je exists, so n x H is continuous as well [67, 110]. With the preceding assump-tions, the interface conditions (2.28) imply that the tangential components of the fields E and H are continuous across any surface in Q,, as is the normal component of B. The normal component n • D may be discontinuous, as a surface electric charge density ge might be induced by the applied source current Je. 2.4 Second-order PDE Formulations It is often not necessary to determine both E and H, as one of these two vector fields can be eliminated. Eliminating H from (2.12) gives the 36 second-order vector Helmholtz equation curl (/x_ 1curlE) - IUJGE = iujJe. (2.29) The corresponding P D E in the time domain is the hyperbolic evolutionary P D E curl (^curl E) + adTE + edTTE = -dtJe. (2.30) After solving for E in either the frequency domain using (2.29) or the time domain using (2.30), the magnetic field intensity H can be recovered if needed using the appropriate form of Faraday's law (i.e. with (2.12a) or (2.4a) respec-tively) . With the other approach, E can be eliminated from (2.12) in the fre-quency domain, giving curl (a _ 1curlH) - iwfiH = curl ( a - 1 J e) . (2.31) Eliminating E in the time domain, however, is less straightforward; dividing by a in (2.31) formally corresponds to computing the operator (a + edt)~l in the time-domain problem. However, with the quasistatic assumption, the term edTE is ignored wherever it occurs, so the time-domain system is curl ( t r - W l U ) + ndtH = curl (a - 1 J e) . (2.32) Notice that (2.32) assumes a ^ 0 which is not necessarily the case in noncon-ducting media such as air. The equations (2.32) are called the equations of magnetic diffusion (see [66, 68]), as they describe a parabolic (rather than hyperbolic) evolutionary 37 system for the magnetic field intensity. The time-domain P D E system for E with the same quasistatic assumption is obtained by dropping the term edttE in (2.30). The frequency-domain systems for E and H under the quasistatic assumption are found by substituting a <— a in (2.29) and (2.31) respectively. 2.5 Weak Formulations Since the classical or strong forms of Maxwell's equations (2.4) or (2.12) implicitly assume high degrees of smoothness in all the fields involved, weak equations posed over suitable functional spaces are required. It suffices to consider spaces consisting of vector fields of finite energy; that is, application of the differential operators grad, cur l , and div should not result in fields with unbounded energy in the standard £ 2 - n o r m . With the domain Q C M.3 given, define the mappings (f,9)a := f f9*dV, (2.33a) (*,i7)n := / t-rfdV, (2.33b) Jn (f,9)an ••= [ f9*dS, (2.33c) J an &»?>«i : = / t-V'dS, (2.33d) J an where * denotes the usual complex conjugate taken pointwise. Suitable norms on spaces of complex scalar and vector fields can be defined in terms of the 38 notation of (2.33a-b), namely, n := v U / ) n , (2-34a) := y/itlh- (2.34b) With the inner products and norms in (2.34) and (2.33a-b), the spaces £ 2 ( f l ) := {/: | | / | |n<oo}, and (2.35a) C2{9) := [£ 2 ( f l ) ] 3 = £ 2 ( f l ) x £ 2 ( f l ) x £ 2 ( f l ) , (2.35b) are Hilbert spaces of complex-valued fields [16, 46]. For variational problems in electromagnetics, the Hilbert spaces £ 2 ( f l ) and £ 2 ( f l ) give rise to the spaces 7i(grad;fl) := { / G £ 2 ( f i ) : g r a d / G C2(Q)} , (2.35c) ?i(curl; fl) := {r? G £ 2 ( f l ) : curlr) G £ 2 ( A ) } , (2.35d) « ( d i v ; fl) := {v G £ 2 ( f l ) : d i v « G £ 2 ( f l ) } . (2.35e) The operators grad, cur l , and divas they occur in (2.35c-e) must be in-terpreted in a weak sense (i.e. as continuous extensions of strong differential operators as they apply to smooth functions; see [16, Chap. 5] for details.) The complex spaces (2.35c-e) are also Hilbert spaces, but not when endowed with the standard £ 2 - inner-products . The appropriate inner product for each 39 space can be inferred from the corresponding Hilbert-space norms: l l / l l ^ g r a d ^ ^ l l / l l n + llgrad/H 2,, (2.36a) ItollTtouruo) : = H"Hn + Wcurlrl\\l ^ d (2.36b) I M I ^ ( d i v ; f 2 ) := H*lln + Udi™||rV (2.36c) All of the spaces (2.35) are Hilbert spaces with their respective norms and inner products in (2.34) and (2.36). Detailed discussions of the required background in functional analysis can be found in [16, 46, 78, 89]). In terms of the above notation, the strong P D E problem (2.23) leads to at least two weak formulations that correspond to equivalent second-order P D E formulations. To see this, given the data e, h G £2(0), define the spaces E e := {£ G ?i(curl; Cl) : n x £ = e on dClE} , (2.37a) Mh := {n G ?i(curl; Cl) : n x 77 = h on dClH} , (2.37b) and their homogeneous counterparts E° := {£ G « (cur l ; Cl) : n x £ = 0 on <9ftB} , (2.37c) M ° := {77 G W(curl ;fi):nxj) = Oon 30^} . (2.37d) The tangential traces n x e and nxh are understood in a distributional sense (£2 fields do not have tangential traces in the strictest sense [46]). With this terminology, the second-order boundary-value P D E curl ((icon)'1 curl IS) - aE = J e + curl ((uop)-1 Jm) (xeCl), (2.38) n x E = n x e (x G dClE), n x {wop)~l [curl-E - Jm] = n x h (x G d f i H ) , 40 leads to the weak problem Find E G W such that, for every £ G E°, ((zu^) _ 1 curl E, cur l£ ) n - (oE, £ ) n = (2.39) (J e, On + J m , cur l£ ) n - (nxh, £)dn . That is, any strong solution E of (2.38) is also a weak solution of (2.39). Notice that the Dirichlet boundary condition for H in (2.23) is replaced by a Neumann boundary condition for E that includes the source field Jm when H is eliminated. Assuming that the frequency u> is not an eigenvalue of the associated homogeneous problem ( ( a ^ - W l £ , c u r l O n - (oE,£)n = 0, the weak problem (2.39) is well-posed [15]. Alternatively, rather than eliminating the magnetic field H, eliminating E gives the boundary-value P D E problem curl ( c r - 1 curl H) — luipH = Jm + curl (o~lJe), (x G fl) V ; V 1 (2.40) n x H = n x h, (x E dflH) n x d~l [curl H — Je] = rt x e, (a; 6 dflE). and the corresponding weak problem Find H G M.h such that, for every rj G H° ( a - 1 curl iJ , curl 77) Q -iuj{pH,rj)n = {Jm,v)n (2-41) + ( a - 1 J e , curl 77) n - (n x e, r))dn . 41 In eliminating either E or H from the original first-order Maxwell's system (2.23), one of the source fields Jm or Je gets incorporated into the boundary conditions and both are combined in the right-hand side of the second-order double-curl PDE. Either of (2.39) or (2.41) can be discretized to find a weak solution describing the desired electromagnetic fields. As is typically for many P D E problems, a single P D E admits numerous weak formulations that are equiva-lent in the event that the solution is sufficiently smooth that a classical solution exists. In a mixed variational formulation, the constitutive relations are not explicitly included in the P D E formulation [20, 21]; this is the approach taken for eddy current and magnetostatic problems in [88, 87]. Weak formulations of boundary-value PDEs are essential for analyses leading to existence of solutions (e.g. [4, 97]). More than being a purely theo-retical tool, weak solutions exist for a number of physically-motivated problems where classical solutions do not. Weak formulations are natural for developing finite-element discretizations that are typically more flexible for dealing with difficult geometries [11, 63, 88]. For the class of geophysical problems studied in this thesis, however, classical P D E formulations leading to finite-volume methods—that are a middle ground between the simplicity of finite-difference methods and the practical strength of finite element methods—are sufficient. 42 Chapter 3 Vector and Scalar Potential Formulations A problem of great importance in geophysics is to compute global elec-tromagnetic fields in three dimensions given a known external source current and an electrical conductivity structure [111]. This is the forward-modeling problem that needs to be solved efficiently within each iteration of a remote-sensing inverse problem [40, 55, 86]. After a thorough description of the forward-modeling problem in Section 3.1, equivalent alternative P D E formu-lations are derived using a Helmholtz decomposition of the electric field. With a new P D E expressed in terms of scalar and vector potential field variables, boundary and interface conditions are presented in Section 3.2, followed by a weak formulation in Section 3.3. 43 3.1 The Forward-modeling Problem The geophysical forward-modeling problem of interest is curl E — iujpH = J , m (ce G fl), (3.1a) curl H - aE = J e (x G fl) (3.1b) n x H = 0 (a; G <9A) (3.1c) with a as given in (2.13). Typically, J m = 0, but the inclusion of Jm in (3.1a) allows for the use of secondary fields (see Section 2.2.6). For most geophysical applications, the dielectric permittivity e is effectively constant, so e = eo, where eo := 8.85 x I O - 1 2 F / m is the permittivity of vacuum. The conductivity a generally varies significantly between the air (a ~ 0 S/m) and the ground (a G [1CT3,102] S/m). The magnetic permeability p may vary over a range of one or two orders of magnitude from po := 4-ir x I O - 7 H/m, the permeability of vacuum. To keep the number of parameters manageable for purposes of the inverse problem, it is common to assume that a is isotropic and piecewise constant. In magnetic ore bodies, assume p is isotropic and piecewise constant like a. For most situations, it is reasonable to assume p = p0 [84, 111]. Finally, in addition to the P D E and boundary condition in (3.1), the interface conditions (2.28) apply at interfaces between distinct material media. The general situation is illustrated in Figure 3.1; the length scale is typically from a few hundred meters to a kilometer. Although an absorbing boundary condition such as (2.25) or even a higher-order absorbing boundary condition may be more physically reasonable 44 Figure 3.1: 2D cross-section of a typical geophysical scenario: e is constant, a and p are piecewise constant. than the homogeneous magnetic (PMC) boundary condition (3.1c), for low-to moderate-frequency ranges (say 0-105 Hz) and for source fields Je and Jm that have compact support in f2, (3.1c) suffices. Under those assumptions, the induced electromagnetic fields decay as |x | —» oo and the dominant behavior of Maxwell's equations is parabolic rather than hyperbolic. This will be the case for an electric or magnetic dipole source, as is common in many geophysical surveys [110, 111]. Thus, for a domain Q that is sufficiently large, assume that a very small fraction of the total energy of the electromagnetic field is stored in the region near the boundary dfl; as such, a homogeneous boundary condition such as (3.1c) suffices for geophysical models. Even when J m and Je have unbounded support but most of the energy of each field (as measured in the jC2-norm) is stored away from dQ, such a homogeneous boundary condition 45 may still be reasonable. The P D E (3.1a-b) can be reformulated in a second-order form by elim-inating the magnetic field H (cf. (2.29)), namely (ztu^curl (pr 1curl.E) - aE = Je + (2w)_1curl (^ _ 1Jm) . (3.2) It is alternatively possible to eliminate E and derive a second-order P D E in the unknown H (cf. (2.31)). For the ranges of conductivities considered here, (3.2) is more favorable; a is much more strongly spatially inhomogeneous than p, the former varying over seven or eight orders of magnitude as compared with two or three for the latter. The dominant term curl (^_1curl E) in (3.2) is far more preferable to the second-order term curl (a~1 curl H) in (2.31) for this reason [53]. In the limit as \ua\ —• 0, the second-order P D E (3.2) reduces to one of the form curl (/z - 1curl.E) = iuiJe + curl (p~l J m ) . (3.3) The curl operator has a nontrivial kernel consisting of any vector field that is a gradient of a scalar field [16, 46]. As such, the P D E (3.3) (with the boundary condition (3.1c)) is ill-posed in the sense that there are infinitely many solutions. Another equation prescribing div E inside Q and a suitable boundary condition is required in order to resolve this ill-posedness; usually, this is the divergence equation (2.12c) with some assumption made about the charge density p [6, 29]. For low frequencies, the displacement current UI€QE is usually very 46 small. Given that the P D E problem (3.1) is ill-posed wherever a = 0, prac-titioners often introduce a small, artificial conductivity, say, a « 10~ 6S/m in those regions [84]. This regularizes the forward-modeling problem (3.1). How-ever, any numerical discretization that faithfully models the continuous prob-lem (3.1) produces an ill-conditioned system of linear algebraic equations, even with the artificial conductivity 0 < a <C 1 in the air. The essential reason is that the small artificial conductivity term can regularize a large, second-order differential term only if it is added to its kernel; this is not the case in (3.1). It is the near-singular nature of the underlying continuous P D E problem that accounts for the ill-conditioned linear systems in experiments (see [84]). Since geophysical data inversions require solving this complex P D E dozens of times, alternative formulations of the analytic problem (3.1) are studied to develop efficient solvers and make computation times feasible. 3.1.1 Helmholtz Decomposition For the time being, questions of differentiability are avoided. The vari-ous electromagnetic fields are assumed sufficiently smooth to allow application of the necessary differential operators. Discussion of boundary conditions and corresponding weak formulations are deferred to Section 3.2 and Section 3.3 respectively. The strategy adopted splits the electric field E into two parts. The first part lies in the range of the curl operator and the second part lies in the kernel of the curl operator. Introduction of the Helmholtz decomposition of E 47 gives E — A + grad (p (3-4) for some vector field A and some scalar field <p [6, 13, 33, 46, 72]. Since curl (grad </>) = 0, the field grad 0 lies in the kernel of the opera-tor curl in (3.1a). With the substitution (3.4) and the preceding observation, when E and H are eliminated, the P D E (3.1a) yields the equation (2a;)-1 curl (^x_1curl A) — d(A + grad</>) = J e The three unknown scalar fields in (3.2) (the components of E) are replaced by four unknown scalar fields (the components of A and the field 4>). Thus, an extra scalar equation is required in addition to (3.5); a useful choice here is the Coulomb gauge condition The field A is typically called a vector potential for the magnetic flux density B since, in the event that J m = 0, (2.12) and (3.4) imply that curl A = curl E = IOJB = tupH. Similarly, the field 4> is called a scalar potential. Thus, the decomposition (3.4),(3.6) yields the new P D E (2u;)-1curl ( / L i - 1 c u r l A ) — o(A + grad 0) = + (iu>) 1 cur l (yU lJ i (3-5) div A = 0. (3-6) J e + (VJ0) lCUX\ (fj, 1Jrn) (3.7a) div A = 0. (3.7b) 48 When ii = fj,0 = constant, the dominant differential term curl ( /Lt _ 1 curl A) in (3.7a) simplifies; in particular, denoting the vector Laplacian (the extension of the Laplacian operator for vector fields) by A , curl ( / /^cur lA) = A o^1 (grad (div A) — A A ) E E - / ^ A A , (3.8) as a result of the Coulomb gauge condition (3.6) [6, 54]. The resulting P D E is [yjjpof1 A A + a(A + grad </>) = — J e — (zw/xo^curl (Jm), (3.9a) d i v A = 0. (3.9b) 3.1.2 Stabilization Generally, the permeability fi is inhomogeneous, so the dominant op-erator curl( /x _ 1 curl) in (3.7a) cannot be replaced with the simpler vector Laplacian operator as in (3.9a). Moreover, the kernel of the dominant oper-ator is nontrivial. A common strategy is to incorporate the gauge condition (3.6) into the equation (3.7a) (see, e.g., [29]). That is, rather than curl (/x _ 1curl (•)), write the dominant operator in (3.7a) as curl (^ _ 1 cur l (•)) — a grad (B div (•)) 49 for some scalar fields a and /?. Given that curl (^ o xcurl (•)) - grad OQ Miv (•)) = 1 A within a medium of homogeneous permeability, the appropriate vector Lapla-cian is recovered when p <— po with the choice a = 1 and (3 = pT1. Hence, the introduction of the vanishing term grad (^ _ 1div A) into (3.7) effectively stabilizes (3.7a); the corresponding system is (iuj)~l (curl (^_1curl A) — grad (/x_1div A)) — a(A + grad^) = The additional term is identically zero for the exact solution of the P D E (3.10), so introduction of the stabilizing term does not change the analytic solution of the problem. This is in contrast with penalty methods where the solution of the penalized problem fails to satisfy some of the constraints [21, 51]. In such cases, selection of a penalty parameter that is large enough to improve conditioning of the discrete equations while small enough to prevent excessive violation of the constraints is a significant challenge (e.g. [51]). 3.1.3 "Pressure-Poisson" Formulation Again deferring the question of boundary conditions, there is a strik-ing resemblance between the constant-yu electromagnetic equations (3.9), the J e + (IUJ) 1curl (/x 1J i (3.10a) div A = 0. (3.10b) 50 Stokes equations Au + gradp = s, (3.11a) div u = 0, (3.11b) and the Oseen equations At* + (w • grad )it + gradp = s, (3.12a) div it = 0, (3.12b) where (3.11) and (3.12) describe steady incompressible fluid flow [46, 51]. In particular, the vector potential field A corresponds to a velocity field u, the scalar potential (b corresponds to a pressure field p, and the Coulomb gauge condition (3.6) corresponds to the incompressibility condition divu = 0 in both (3.11) and (3.12). The P D E (3.9) for electromagnetics and the corre-sponding PDEs (3.11) and (3.12) for fluids differ in the (complex) Helmholtz shift IUSOA in (3.9a) and the advective term (w • grad)w in (3.12a); bar-ring these differences, the systems (3.9), (3.11), and (3.12) have very similar structure in the dominant differential terms for (A, <p) and (u,p) respectively. Motivated by the analogy with fluid flow, take the divergence of (3.1b) and replace the Coulomb gauge in (3.10) with the resulting equation. The P D E obtained is (zo;)-1 (curl (yu-1curl A) — grad ( p : - 1 d i v A) ) — a(A + grad</>) = J E + (iu>) xcurl (fi 1Jm) (3.13a) div (<T(A + grad </>)) = —div J E . (3.13b) 51 The equation (3.13b) resembles the pressure-Poisson equation (PPE) used in incompressible flow [51, 99, 108]. The P D E (3.13) is used in [13] to derive a finite-element discretization of a PPE-type systems in (A, <p) for a large scale geomagnetic problem where p = po. The same P D E is studied in [6, 54] to yield finite-volume methods for other low-frequency geophysical applications. In [53], a discretization of a related PPE-type system in (A, </>) is derived for the more general case of inhomogeneous permeability with the additional stabilizing term. Solving direct discretizations of the indefinite P D E system (3.7) is pos-sible by techniques analogous to those in [38, 39, 51] for the fluid-flow problems (3.11) and (3.12). The principle distinction between (3.7) or (3.10) and the fluid-flow problems (3.11) or (3.12) is that the decomposition (3.4),(3.6) is ar-bitrary (although A and <p do have physical interpretations [6, 13]). There is, to a degree, more freedom in prescribing boundary conditions for the sys-tem (3.13) for this family of geophysical applications. This freedom contrasts the additional problems with boundary conditions arising in pressure-Poisson formulations for Stokes' flow [51, 99]. The left-hand side of the system (3.13) includes the linear differential operator If ui is not too large in magnitude, the differential operator in (3.14a) is di-(iu) 1 Ap,-a I -cr grad (•) (3.14a) A ^ := curl (p ^url^)) — grad (/n 1div(-)) . (3.14b) 52 agonally dominant; that is, the leading order differential terms occur on the main diagonal. As a result, the corresponding scalar PDEs constituting (3.13) are weakly coupled. When (3.13) is discretized, the diagonal dominance of the operator leads to matrices with block-diagonal dominance; this property is used to advantage in devising iterative solvers. By contrast, standard for-mulations of Maxwell's equations (such as (3.2) for E) have strong coupling between the scalar components. The poor convergence behavior of iterative methods applied to the corresponding discretizations of (3.2) relates to this coupling (see [53, 84]). 3.2 Boundary Conditions To close the stabilized forward-modeling P D E (3.13), suitable boundary conditions for the unknowns (A, cj)) must be imposed. The key steps taken in the preceding derivation are: • introduction of a Helmholtz decomposition with a Coulomb gauge; • introduction of a stabilizing term; • differentiation to obtain a "PPE" to replace the gauge condition. The original P D E (3.1) has a single boundary condition and the final boundary-value problem has two additional boundary conditions derived below. It is convenient to use auxiliary variables to derive the boundary con-ditions. Hence, introduce the generalized current density j:=-ZE (3.15) 53 into (3.1). In terms of j , (3.1b) becomes curlH + j = Je. (3.16) Eliminating H from (3.1) and applying the definition (3.15), the boundary-value problem (3.1) becomes (2w) - 1curl ( / / _ 1 cur l£) + j = F, (3.17a) j + aE = 0, (3.17b) together with the magnetic boundary condition n x / j " 1 c u r l £ = n x / i . (3.17c) For convenience, the source and boundary terms in (3.1) are written as h := (p-1 Jm), (3.17d) F := Je + [VJJ)-1 curlfc. (3.17e) The transformation (3.15) is taken pointwise and involves no derivatives, so the boundary conditions remain unaffected in (3.17). The substitution of the Helmholtz decomposition (3.4),(3.6) into (3.17) requires another boundary condition. It is well known that choosing / Jo. n • A = 0 on dQ and (3.18a) <j>dV = 0 (3.18b) is sufficient to guarantee that the Helmholtz decomposition is unique [46]. It is also possible to partition the boundary, prescribing (3.18a) on one portion of 54 dfl and prescribing a Neumann-type condition on <p over the remainder of dfl [90].' In either case, clS with the Neumann problem for Poisson's equation, the scalar field 4> is determined up to a scalar constant, so either a point condition or an integral condition like (3.18b) is required to fix the field <b [30, 89, 103]. The new boundary-value problem, then, is the P D E (ttj)_1curl 0 _ 1curlA) +j = F, (3.19a) d i v A = 0, (3.19b) j + a(A + grad </>) = 0, (3.19c) with boundary conditions and constraint n x (/^curl A) = n x h, (3.19d) n-A = 0, (3.19e) [ (bdV = 0. (3.19f) Since the B V P (3.17) is known to be well-posed [15] and since the Helmholtz decomposition with (3.19e,f) is also known to be unique, it follows that the B V P (3.19) also has a unique solution. The dominant operator in the equation (3.19) is stabilized through sub-traction of a term that lies in the kernel of the curl operator, namely the term grad ( / i _ 1 div A). Rather than expressing the stabilized P D E as in (3.10a), introduce a new variable ^ : = ^ _ 1 d i v A . (3.20) The stabilizing term in (3.10a), then, is grad^, and the equation (3.19a) 55 becomes (UJ)-1 (curl (/Li_1curl A) - grad ip) + j = F. (3.21) If the Coulomb gauge condition applies, ip = 0 and the stabilizing term grad ip in (3.21) also vanishes identically. Next, the divergence is taken of the equation (3.19a), yielding divj = div Je. (3.22) Replacing the Coulomb gauge condition with (3.22) gives (lu)'1 (curl (^ -1curl A) - grad rp) +j = F, (3.23a) divj = div J e , (3.23b) j + a{A +grad <j)) = 0, (3.23c) pip - d i v A - 0 , (3.23d) with the two boundary conditions and the constraint (3.19d-f). The extra differentiation to get (3.23b) requires an additional boundary condition to close the P D E problem (3.23). To derive the appropriate additional BC, write Ampere's law (3.1b) as curlH + j = Je (3.24) Assuming (3.24) holds on a small surface S C dfl, then, integrating the normal projection of (3.24) onto S (with the usual outward orientation) gives / H-rd£+ f n-jdS= [ n-JedS. Jas Js Js 56 By (3.1c), the first integral vanishes, leaving Js Js [ n-jdS= [ n • JedS (3.25) for any surface S C dfl. It follows that n • j = n • Je on dfl. (3.26) The procedure above of integrating the normal component of Ampere's law corresponds exactly to the derivation of the pressure boundary condition for the pressure-Poisson equation as outlined by Gresho and Sani [50]. The importance of the boundary condition (3.26) is that, as in the con-text of incompressible fluid flow, it is consistent with the requirement div A = 0 on dfl [51]. Take the divergence of (3.23a); by (3.22), that is, ip satisfies Laplace's equation. Writing the stabilized equation (3.21) as and applying the same arguments used in deriving (3.26), it follows that ip satisfies a homogeneous Neumann boundary condition, i.e., dip/dn = 0 on dfl. Hence, ip is spatially constant; an additional constraint such as Jo. pins down ip uniquely. The Coulomb gauge condition (3.6) is thus recovered throughout fl. Atp = 0; (3.27) curl H — {iuj) :grad ip + j = F (3.28) (3.29) 57 Hence, the final boundary-value problem is (iu)-1 (curl ( / / _ 1 curl A) - grad ip) + j = F, (3.30a) divj = div J e , (3.30b) j + a(A + grad <p) - 0, (3.30c) pip - d i v A = 0, (3.30d) with boundary conditions n x (/u - 1curl A) — n x h, (3.30e) n- A — 0, (3.30f) n • j = n • Je. (3.30g) and global constraints / ipdV = 0, (3.30h) [ (pdV = 0. (3.30i) An equivalent boundary-value problem in terms of the variables A and 4> alone is the P D E (IOO)'1 (curl (^ _ 1 cur lA) — grad ( / i _ 1 divA)) — a (A + grad0) = J e + (zw) _ 1curl (/x_1 J m ) , (3.31a) div [a(A + grad 0)] = -d iv J e , (3.31b) 58 with boundary conditions n x (fi x cur l A) = n x (/i 1JTn) (3.32a) n • A = 0, (3.32b) n • CT(A + grad <p) = —n • Je. (3.32c) and constraints / fi^divAdV = 0 (3.32d) f <f>dV = 0. (3.32e) The boundary conditions (3.32a,c) involve the material coefficients; in the event that the material coefficients are in fact tensor fields, the interactions of cr, e, and fi with the normal field on dfl are still accounted for by (3.32). For the simple case of isotropic tensor fields considered here, some further simplifications can be made. In particular, n x curl A = n x Jm, (3.33a) For most geophysical applications, it is undesirable to have action of the source fields at the boundary. If the source fields have bounded support, then the boundary conditions (3.33) are homogeneous. 3.2.1 Interface Conditions With the introduction of different electromagnetic field variables, the new fields satisfy different jump conditions at the interfaces of conducting (3.33b) 59 media. The same procedure used to derive (2.28) (cf. Section 2.3.2) yields the conditions n i r „ i [A1 - A11] = 0, (3.34a) n 1 x [A 1 - An] = j m , (3.34b) n 1 • [e\Al + grad cp1) - e n ( A n + grad 011)] = gei (3.34c) [j1 - ju] = n1 • [Jel - J*] (3.34d) n1 ^ across any surface T c A (cf. Section 2.3.2). The surface magnetic current density Jm is often identically zero; where Jm ^ 0, it generally does not have surface layers. As a result, n 1 x [A 1 - Au] = 0, (3.34e) i.e., the tangential component of A is continuous across T. Similarly, with the assumption e = eo in fl, (3.34a,c) imply nl • [grad cp1 - grad cp11] = ge/e0. (3.34f) Finally, provided the surface T does not coincide with a surface of discontinuity of the applied electric source current J e , the right-hand side of (3.34d) vanishes also, i.e., n1 r - i [j1-jll]=0. (3.34g) The preceding assumptions are not particularly restrictive in the context of this work or even in more general contexts. The resulting interface conditions (3.34a,e-g) are of great utility in deriving discretizations in Section 4.3. 60 3.3 Weak Formulations There are naturally many weak problems equivalent to the P D E prob-lem (3.31) that can be examined in a variational context. However, beyond pointing the way to different numerical schemes and providing a basis for finite-element analysis, weak formulations can provide insight into the interplay of boundary conditions with the differential operators. Although demonstrating the existence of solutions from weak formulations falls outside the scope of this work, a brief discussion of one particular weak formulation (similar to that applied in [13]) is included for completeness. In terms of the notation (2.33) introduced in Chapter 2, the divergence theorem implies (curl £, ri)n - (£, curl rj)a = (nx£,r])dQ, (3.35a) (u, div V)n + (grad u,V)a = (u,n-V)dn (3.35b) for any £, r] G 7i(curl ;fl),V e ?i(div ; fl), and u G 7Y(grad ; fl). Multiplying (3.31a) by some test vector field £ and integrating by parts with (3.35) gives the weak equality (iu)-la^(A, i) - (a(A + g r a d 0 ) , £ ) n - {iTldxv A,n • £)QQ = f(£) + (zu)-1(nx(p-1(Jrn-curlA)),Z)dn, (3.36a) where a^A,£) := (/z -1curl A, curl £) n + (^ _ 1div A , d i v £ ) n , (3.36b) / ( £ ) := ( J e , O n + M " 1 ( / i - 1 J m , c u r U ) n . (3.36c) 61 Similarly, multiplying (3.31b) by a scalar test function v and integrating gives (a{A + grad <f>), gradv)n = - (J e, gradw)n + (n • (o(A + grad0) + Je), v)9n . (3.36d) For a solution of the P D E (3.31), the boundary integrals on the right-hand sides of (3.36a,d) vanish by virtue of (3.32a,c); that is, the boundary conditions (3.32a,c) are natural boundary conditions for this choice of variational form. The term (/x_1div A,n • £)an on the left-hand side of (3.36a) does not vanish unless some constraint is imposed on the space of functions from which the trial and testing functions are drawn. Assuming n • £ = 0 for all admissible vector fields £, then the boundary integral term in question vanishes. Thus, imposing (3.32b) is an essential boundary condition. Combining the above arguments, define the spaces A 0 := {£ e « ( d i v ; fl) D ?i(curl; fl) : n • £ = 0 on dfl} , (3.37a) $ := 7i(grad;fl). (3.37b) A variational problem that is equivalent to the B V P (3.31) subject to boundary conditions (3.32) is to find (A, (j>) G A 0 x $ such that (iu)-X(A,Z) - (a(A +grade/)),On = / ( 0 V£ G A 0 , (3.38a) (d{A + grad 0), grad v)n = - {Je, grad v)n Vv G (3.38b) A strong solution of the classical BVP (3.31),(3.32) will also satisfy the varia-tional problem (3.38). Inserting a test function £ = grad?; G A 0 into (3.38a) and using (3.38b), it is clear that a solution of the weak variational problem is also weakly divergence-free. 62 Chapter 4 Finite-Volume Discretizations Having chosen an analytic formulation ((3.31) with (3.32)) for Maxwell's equations, the next step is to select a suitable discretization. While certain applications call for finite element discretizations in order to capture the ge-ometry accurately (see [16, 63, 88, 87] among many others), for the family of geophysical applications of relevance here, finite-volume discretizations suffice. After presentation of a description of the discrete domain in Cartesian coor-dinates, the well-known Yee scheme is derived using standard finite-volume arguments. A finite-volume discretization of the P D E (3.31) is subsequently presented in Section 4.3; this discretization extends the ideas of the Yee scheme to account for inhomogeneous media in a logical manner. Finally, Section 4.4 shows some analysis and experiments that justify the choice of coarse discrete representations of the inhomogeneous media. 4.1 Domain of Discretization Assume the domain of the continuous problem (3.31) in Cartesian co-63 ordinates is the rectangular box (0, Lx) x (0, Ly) x (0, Lz) C K 3 . The various finite-volume discretizations are derived using a tensor-product grid of Carte-sian coordinates whose coordinates correspond to points in space at which the discrete grid functions approximate the appropriate continuous functions. The same grid and notation in this section underlies the Yee discretization in Section 4.2 and the discretization of the desired system (3.31) presented in Section 4.3. Consider first the grid in the x-direction. There are Nx subintervals partitioning the interval (0, Lx). The widths of the subintervals are denoted hx (i = 1,..., Nx) and the centers of the subintervals are xt (i — 1,..., Nx). The end-points of the subintervals have half-integer indices, namely x 4_ 1 / 2 (i = 1,..., Nx + 1). The dual grid comprises the subintervals with end-points xt of widths h*_i given by hi i : = X i - x^ (i = l,...,Nx + l). (4.1) The ghost points x0 := xt - h* and xNx+l := xNx + hxNx (4.2) may be included to facilitate derivation of boundary conditions; any discrete unknowns associated with grid points are eliminated using discrete approxi-mations of the boundary conditions (e.g. [44, 99, 113]). Assume that the intervals (0, Ly) and (0, Lz) are partitioned into Ny and Nz subintervals respectively with dual grids defined analogously. The vertices (xi±%> Uj±^i zk±\) are the corners of the cells of the primary grid; that is, the 64 boxes Vi,j,k •= (x^i,xi+i) x x (zk_i,zk+i) (4.3a) (i = l,...Nx-j = l,...Ny-k = l,...Nz) with faces S?±li,k := ixi±0 x (Vj-l'Vi+l) x (4.3b) ^ i i . f c : = ( ^ - ^ + 1 ) x { % ± ^ } x (zk-$,Zk+i), (4.3c) ^ , * ± i : = fo-i>*«+i) x x ( 4- 3 d) (z = l , . . . / V x ; j = l , . . . iV 2 / ;A: = l ) . . . 7V z ) constitute the cells of the primary grid. The dual grid comprises the cells between the centroids of the primary grid cells; that is, the boxes (4.4a) Nz-1) (4.4b) (4.4c) (4.4d) Nz - 1). Notice that, while the centroids of the cells and faces of the primary grid coincide with points in the dual grid, the centroids of the cells and faces of the 65 V i + i J + i i k + i := (xifxi+1) x ) x (zk,zk+1) {i = l1...Nx-VJ = l,...Ny-l-k = l,... with faces S - + i ± i J + i M i := {Zi + i ± l } x (&>&+0 x (^,zk+1), Si+lj+i±ijk+i '•= (xitxi+1) x { % + i ± i } x (zk,zk+1), S-+y+i_tk+i±i_ := {xuxi+1) x (y,,y,+ 1) x {zk+i±i} (i = l , . . . 7 V B - l ; j = l , . . . / V 1 / - l ; f c = l , . . . dual grid do not necessarily coincide with points in the primary grid unless the grid spacings are uniform [80]. The basic geometry of the grids is illustrated for the case of uniform grid spacing in Figure 4.1. Primary grid Dual grid Figure 4.1: 2D cross-section of the primary and dual grids (uniform spacing). The corners of the dual grid are the cell centers of the primary grid and conversely. 4.2 The Yee Discretization Yee [115] derived a set of finite-difference equations approximating Max-well's equations in the time-domain. In this scheme, both the electric and magnetic fields E and H are discretized and are explicitly marched forward in time using the coupled curl equations (2.4a-b) as opposed to discretizing E alone and using the second-order wave equation formulation (2.30). The finite-difference method derived is often called the Yee algorithm or the Finite-66 Difference-Time-Domain algorithm (FD-TD) in the engineering literature (as in [104, 105]). The related spatial discretization in the frequency domain is the Finite-Difference-Frequency-Domain algorithm. In this thesis, the termi-nology "Yee algorithm" is avoided in favor of the more accurate terminology "Yee discretization" or "Yee scheme." Yee's important contribution lies in rec-ognizing how to discretize Maxwell's equations in space in conservative form; from the perspective of evolutionary ODEs or PDEs, the algorithmic aspects of the "Yee algorithm" come from well-established time-stepping techniques (e.g. [8, 91]). The basic idea for the Yee discretization is to center the components of E and H in three-dimensional space so that each component of E is sur-rounded by four circulating components of H and conversely (except at the boundaries). With the unknown fields staggered on a grid, Faraday's Law and Ampere's Law are invoked and approximated along the contours. That is, the circulation of E in a loop is linked with the flux of the magnetic field H through that loop while the circulation of H in a loop on a dual lattice is linked with the displacement current flux of the electric field E through that loop. For the engineering applications for which Yee's discretization was origi-nally derived, there are a number of practical advantages. Firstly, the locations of the components of E and H in the spatial lattice implicitly enforce discrete versions of the divergence conditions (2.4c-d). In fact, although originally de-rived as a finite-difference scheme, the Yee discretization can be derived as a 67 finite-volume method. Thus, the conservation properties arise as a common consequence of finite-volume techniques; for certain application families, it is important to preserve such global conservation laws [66, 104]. Secondly, the resulting finite-difference approximations for the space derivatives used in the curl operators are centered-difference approximations and hence are second-order accurate. Thirdly, when the staggered-grid spatial discretization is cou-pled with a leap-frog discretization in time, the resulting time-stepping method is explicit. The resulting discrete electromagnetic field can be evolved with-out solving large sparse linear systems at each time step. This can be useful for simulations of waveguides or other three-dimensional electronic structures over long time scales, since solving the linear systems involved would be pro-hibitive [104]. Finally, since the discrete fields are specified in terms of normal or tangential components, specification of certain boundary conditions can be simplified. For instance, given the tangential trace n x e o f some known field e, if the boundary condition n x E = n x e on dfl is given, the tangential components of the discrete electric field Eh on the discrete boundary of a Cartesian tensor-product grid are given directly in the discretization. 4.2.1 The Discrete Fields For the Yee discretization, the various components of the vector fields are staggered in space to approximate Maxwell's equations in a way that au-tomatically ensures the integral conservation laws (2.4c-d) hold on a discrete level. In particular, the components of the discrete electric field Eh are defined 68 at the centers of the edges of the cells of the primary grid and are tangential to those edges. Explicitly, if E = (Ex,Ey,Ez)T in Cartesian coordinates, the discrete electric field unknowns are EXi,j±lk±i, - Ex(xi,yj±i,zk±h), (4.5a) EVi±y,k±\ ^Ey(xi±i,yj,zk±i), (4.5b) EZi±lj±ik - Ez(xi±i,yj±i,zk), (4.5c) where i — 1,... Nx, j = 1,... Ny, and k = 1,... Nz. From this, it follows that, for example, dEx E x i i + i k + i - E x i i + i k _ i and similarly for other spatial derivatives of the components of E; that is, at the centers of cell faces, centered finite-difference approximations of derivatives of components of E are second-order accurate. The components of the discrete magnetic field Hh are defined at the centers of the faces of the cells of the primary grid and are normal to those cell faces. Thus, given H = (Hx, Hy, HZ)T, the magnetic field discrete unknowns are HXi±lj,k - Hx(xi±h,Vj,zk), (4.6a) H ' ^ k - H ^ x ^ ^ z , ) , (4.6b) Hzi,j,k±i - Hz(xuyj,zk±i), (4.6c) with i = 1,... Nx, j = 1,... Ny, and k = 1,... iV2. The placement of the various discrete solution components on the cell is depicted in Figure 4.2 and Figure 4.3. 69 Figure 4.2: Staggering of components of the discrete electric field Eh for the Yee scheme. 4.2.2 Derivation of Yee Scheme The approach presented here uses finite-volume arguments along the lines of [27] or [80] as opposed to finite-difference arguments such as in [104] or Yee's original work [115]. Although the Yee scheme was originally derived on a uniform grid, the presentation here allows for nonuniform grids. Also, the frequency-domain Maxwell equations are discretized (as in [77, 84]) rather than the time-domain equations (2.4) used more often in engineering applica-tions. Furthermore, for generality, a magnetic current source J m is included in the right-hand side of the discretization of Faraday's law, as in (2.16a); this generalization may be helpful for magnetotelluric models as in Section 2.2.7. Finally, for the moment, assume that the material parameters cr, e, and // 70 Figure 4.3: Staggering of components of the discrete magnetic field Hh for the Yee scheme. are all homogeneous and isotropic, i.e. scalar constants. This assumption is relaxed in Section 4.3 where the discretization of the system (3.31) is derived in terms of discrete vector and scalar potentials in inhomogeneous, isotropic media. To approximate Maxwell's equations (2.16) on the nonuniform grid flh, consider integrating the Maxwell curl equations (2.16a-b) as in [80, 104]. For any smooth surface S C fl with boundary dS, Stokes' theorem and Faraday's Law (2.16a) give I E-rde-iuj [ nH-ndS= [ Jm-ndS. (4.7) Jes Js Js The various discrete equations that correspond to the components of (2.16a) 71 are derived by choosing S as a face of the primary grid. First, take S = Sf_idik := {x^i} x (yj_i,yj+i) x (zk_i,zk+i) in (4.7) as given in (4.3b). Approximating the line integral with four quadra-ture points at the midpoints at the center of each line segment and approxi-mating the surface integrals with a single quadrature point at the centroid of S gives - ^ j f W . ^ = hyKFmi_^k. (4.8a) The remaining two discrete equations are obtained using from (4.3c) and S = SfJjk_i := ( x ^ i . x ^ i ) x (yj_i,yj+i) x {zk_i} from (4.3d) in (4.7); the corresponding finite-volume equations for the mag-netic field Hh are - uj^hiHy^j, = KKjymh3_hk, (4.8b) - l u n h W H ' ^ = Khyym^k_h. (4.8c) 72 To discretize the electric field equations, use the integral form of (2.16b), namely I H • T d£ — f aE-ndS= [ Je-ndS. (4.8d) Jas Js Js The discrete equations for the electric field arise from choosing the surfaces 5 as the faces of the dual grid formed by joining the centroids of the cells of the primary grid. For instance, choosing S = S.._ik_i := {xt} x iVi-uVi) x (zk_uzk) from (4.4b) in (4.8d), and using one-point quadrature formulae to approximate the integrals gives - h]_hK_hdBr^M_h = / » ; . 4 ^ _ 4 J T y _ i l f c _ 4 . (4.8e) For the corresponding equations in the y- and z-directions, substitute 5 = 5 ? , . . i : = ( X i _ i , a ; i ) x {yj} x (zk_uzk) as in (4.4c), and l ~ 2'J<K 2 5 = 5f i . _ i k := ( i ^ . X i ) x {y^uVj) x { zJ as in (4.4d) ' 2 " 2' into (4.8d), and repeat the procedure outlined above, respectively giving - h^hUfW^y^i = hUhl_hJyei_Uk_^ and (4.8f) - /if rfE'^i , i fc = /if i/iw x J eVi 7-_i fc- (4-8g) 73 For the cells adjacent to the boundary dfl, integrating (2.16a) with S C dfl gives a discrete equation that can be simplified by appropriate boundary con-ditions (see [99] for examples). With the discrete fields situated as they are on the staggered grids, the Yee scheme is second-order accurate on a uniform grid as can be determined by simple Taylor expansions [91, 104]; for greater accuracy on nonuniform grids, typically the sampling of the source terms has to be done more carefully (e.g. [7]). An alternative convergence analysis by Monk and Suli establishes the second-order accuracy of the discretization on nonuniform grids in suit-able discrete energy norms [80]. The above finite-volume approach can also be extended to include systems involving inhomogeneous media and nonorthog-onal or unstructured grids [66, 104]. Finally, highly accurate finite-volume discretizations of Maxwell's equations can also be derived from conservation-law form (e.g. [31]); such discretizations become more important in hyperbolic high-frequency transmission problems where accuracy is crucial. 4.3 Discretization Using Potentials The goal now is to introduce a finite-volume discretization of the P D E system (3.31) derived in Chapter 3. For the case of p = po = constant, a finite-volume discretization is derived in [54]; the presentation below considers the case where /x varies and is based on the original derivation in [53]. Piecewise continuity of the isotropic, inhomogeneous material tensor fields is assumed up to the resolution of the spatial grid. 74 The domain 0 is subdivided into cells as outlined already in Section 4.1. Within each cell Vijtk of the primary grid, the conductivity a, the permeability fi, and the permittivity e are all assumed constant; that is, 0" = crj^ fc, n = A*j,j,fc, and e = e i j i f c (4.9) throughout the cell Vijtk- However, a, //, and e may be discontinuous between adjacent cells. The material fields could be varying smoothly within cells, but it is simpler to assume that the material properties are constant in each cell. The domain, then, is composed of blocks of materials of distinct electric and magnetic properties. In practice, the variations in e are much smaller in magnitude than those in a and p for geophysical problems (see Section 3.1). For the purpose of deriving a discretization of (3.31) in terms of discrete grid functions (Ah,<fih), it is advantageous to use the auxiliary or intermedi-ate variables used in Section 3.2. Using the additional variables to derive a discretization is akin to mixed finite element methods that are commonly applied to problems with highly discontinuous coefficients [20, 21, 88]. This approach is employed in both [53] and [54] although the final discrete system is a discretization of (3.31) with discrete approximations (Ah, 4>H) of (A, 0) only. Thus, in terms of the variables introduced in Chapter 3, the system 75 under consideration is curl H - (ico)~lgrad ib + j = Je, (4.10a) d i v j = d i v J e , (4.10b) fii/) - div A = 0, (4.10c) j + o (A + grad cb) = 0, (4. lOd) curl A — lupH = Jm (4.10e) 4.3.1 The Discrete Fields Having decided to use the variables (A,j,H,<p,ip) to derive a finite-volume discretization, it is necessary to decide where on the staggered grid to place the components of the discrete approximations (Ah,jh,Hh,<ph,tph). This element is left out of the derivation of the Yee discretization in Section 4.2. There, the components of Eh and Hh are chosen as in (4.5) and (4.6) to en-sure the discrete divergence conditions hold. It is equally reasonable, then, to reverse the roles of Eh and Hh on the staggered grid, yielding a discretiza-tion with Hh tangential to cell edges and Eh normal to cell faces; within homogeneous media, the choice is arbitrary. With inhomogeneous media, certain field components are known to be continuous or discontinuous across material interfaces. It is the interface con-ditions (2.28) together with the discretization (4.9) of material coefficients that points towards more robust discretizations. Under the assumptions listed in Section 3.2.1, the interface conditions across cell surfaces for the variables A, 76 j , H, and <p are as in (2.28) and (3.34), i.e., n • [A1 - A11} = 0, (4.11a) n A1 - A11 = 0 (4.11b) n •[j1-3U]=0 (4.11c) n x [H1 - Hu] = 0 (4.11d) n • [grad^1 - grad^11] = ge/e0, (4.11e) where ge is a surface charge density. The conditions (4.11) imply that n • j is continuous across material interfaces while n • E is not. Moreover, n • grad (p inherits the discontinuity in n • E, while A is continuous. Motivated by the conditions (4.11), the locations of the discrete un-knowns become clear. Since discrete divergences are computed by integrating over cell volumes, the discrete unknowns Ah and jh corresponding to the vec-tor potential A and the generalized current field j respectively are prescribed on cell faces of the primary grid with components normal to cell faces. By (4.10c), the unknowns comprising the discrete field iph should be located at the centers of the cell volumes since that is where the discrete approximation to div A is centered. Similarly, the discrete field <ph sits at cell centers; with that choice, the approximation of grad (p is normal to the centers of cell faces. Since, from (2.28d), the tangential components of the magnetic field H are continuous across interfaces in the absence of surface currents je, the discrete unknowns Hh are set at the centers of the edges of cells of the pri-mary grid with components tangential to those edges. This contrasts with the 77 specification of Hh in the Yee discretization as given in (4.6); Hh is tangential to cell edges in this scheme while it is normal to cell faces in the Yee scheme. Then, in Cartesian coordinates, the discrete approximations of A = (Ax,Ay,A*)T, j = (f,jy,jz)T, H = (H*,Hy,H*)T, cb, and V are Axi±i)jtk^Ax{xi±h,yj,zk), (4.12a) Ayij±iJ.*Ay(x«V,±h,z*), (4-12b) A\j,k±\ - A z i x u yj,Zk±l), (4.12c) A ± ^ - , f c - J x ( ^ ± i » ^ . ^ ) > ( 4 - 1 2 d ) 3Vij±lk-3v(xi>y*l,z<.)> (4-12e) 3ij,k±±-3'(Xi,yi,zk±h), (4.12f) ^ i ^ f e f e j . ^ ) - (4-l2g) Hyi±i,jtk±^Hy(xi±h,y,,zk±h), (4.12h) H z i ± i 2 , j ± h k ^ H z ( x i ± i , y j ± h , z k ) , (4.12i) 0ij,fc - 0(Zi, ( 4- 1 2J) il>ij,k-il>{Xi,yj,zk), (4.12k) where z = l , . . . i V x , j = l,...Ny, and A; = 1,...NZ. The discrete source field Jg is normal to cell faces as is A while the discrete source field Jhm lies tangential to cell edges like Hh. 4.3.2 Derivation of the Finite-Volume Scheme The strategy to derive the discrete system is as follows: 78 1. Integrate (4.10a) over the faces of primary-grid cells. 2. Integrate (4.10b) over the primary-grid cells. 3. Integrate (4.10c) over the primary-grid cells to get discrete equations for 4. Integrate (4.10d) over lines between centers of adjacent primary-grid cells to get discrete equations for j h. 5. Integrate (4.10e) over the faces of dual-grid cells to get discrete equations for Hh. In the last three steps, the discrete equations derived are used to eliminate the discrete fields iph, j h, and Hh from the discrete equations derived in the first two steps. Where appropriate, the boundary conditions (3.32) can be applied to close the resulting discrete equations [44, 99, 113]. The end result is a large linear system that involves only the unknown discrete fields Ah and <ph. To begin, integrate the stabilized equation (4.10a) over faces of the primary grid as in the derivation of the Yee scheme. For any surface S C A, (4.10a) together with Stokes' Theorem gives / H • T d£ — (IOJ)"1 [ gradi />-ndS+ [ j-ndS= [ JendS. (4.13) Jas Js Js Js The surfaces to use in (4.13) are the faces of the cells of the primary grid (the rectangles Sx x .,, Sy. t ., and Sz., 1 from (4.3)). Substitute S = Sx_, .. into (4.13) and use a centered finite-differences to approximate dxtp. Then, 79 with the same discrete approximations as in (4.8a) for the remaining integral terms, the discrete equation resulting is + m + = VKJ:^. ' (4.14a) The corresponding equations in the y- and ^-directions are obtained using S = £>lj_i k a n d S*jk_i respectively in (4.13); the resulting discrete equations are K (Hxi,j-lk+i ~ Hxitj_ik_i^ - hzk (H'i+ij_itk - Hzi_i + K K = K K j y ^ (4.i4b) Next, for any volume V C fl, (4.10b) gives (p j • ndS — d> Je • ndS. Jav Jav Choose V = Vij,k from (4.3a) to derive h)hzk (fi+ijt - fi-y,k) + KK (fij+^k ~ fij-ik) + K K ^ k H - J : ^ . (4.i5) 80 Having derived (4.14) and (4.15), the remaining equations are integrated to eliminate the discrete unknowns iph, jh, and Hh. The resulting system involves Ah and (bh only; it is this system that is studied in subsequent chapters. Eliminating iph is straightforward; by (4.10c), for any volume 7 c O , / fiipdV - / div Adv = 0. Jv Jv Since /x = Hijyk is constant in cells, as in the derivation of (4.15), choose V = ViJtk, yielding + ^ { A ' i J M h - A ' i J J l _ h ) ] . (4.16) Since the field grade/) is potentially less smooth than <7grad0 across an interface between primary grid cells, it is preferable to express (4.10d) by dividing through by a first as in [54]. Then, for any path F C fl, J [a'1 j + A + grad cp] -Tde = 0. (4.17) Choosing T in (4.17) to be a line between cell centers such as T = %_iJik •= (Xi-i, xt) x { j / J x {zk} , the line integral can be approximated by [ £ j l i £ £ + A \ _ y \ + fa* - b-ij* = 0. (4.18) In (4.18), the harmonic average of a across the adjacent cells V j - i j ^ and Vijtk is defined by _hX (r dx x - 1 0 81 By (4.9), the integral above evaluates to / hx hx \ _ 1 at_,ik = hx ! + . (4.19a) Similar averages taken in the y- and ^-directions give du-l* = h U ( i t + W J - ) ' ( 4 1 9 b ) The derivation of harmonic—rather than arithmetic—averages is natural in this context upon application of numerical considerations (viz. integration of rough quantities). Using harmonic averages becomes important in practice when the coefficient a may jump by several orders of magnitude [54, 96, 113]. The discrete approximation (4.18) can be solved to yield an equation for j X i _ i J t k , namely f , - h j , k = ( A ' H * + 6^~f;-hj^ • . (4.20a) The integration of (4.17) along lines in the y- and ^-directions yields iVi,* = - 9 « - J * + , (4.20b) = ( A ' ^ + 0 ' " \ " / U t - ' ) (4.20c) for j v i ^ _ \ k and j z i i k _ i respectively. Thus, using the definitions (4.19), the equations (4.20) can be used to eliminate jh in (4.14) and (4.15). It remains to eliminate the magnetic field unknowns Hh from (4.14) by integrating the equation (4.10e). Since the tangential components of the field 82 B = pH are discontinuous while the tangential components of the field H are continuous across interfaces (see (2.28)), it makes more sense to integrate pH in (4.10e). Thus, the equation (4.10e) implies iu j pH -ndS = - I Jm-ndS+ I A-rdl (4.21) JS Js JdS for any surface S. The surfaces of integration to use in (4.21) are the faces of the dual grid, namely Sx. i , i, Sy 1 ., j, and Sx i . r , from (4.4b-d). Taking Sx. i , i, use single quadrature points to approximate the integrals above, so K.h - ^Vi,*-*) - fcj-j - ^ Vi*-i) • (4-22) In (4.22), the value of JJL on the surface Sx. i , i is approximated by the 1 '3 2 ' 2 ordinary average ^ j - i ^ - i := (^_i^ Li) ^ p(xi,y,z)dydz. (4.23) Assuming the permeability is piecewise constant, this simplifies to = £ ( 4 - 2 ^ 9 , r=±± 83 with corresponding averages / V l j . f c - 1 = ( 4 /C^Ll) 1 E hU+,hl-h+rtii-$+PJ,k-h+r> (4-24b) ^•-i^Kl^l)"1 E ^-i^ J-i^ -i^ -i^ (4.24c) in the and z-directions respectively. The discrete equations for H are = (VKj-ik-i)-1 (A'id^~A'iJ-1^ (4.25a) Ayi _ i fc - Ayi _ i fc.j \ - " 2 ' \ , t ' " - > (4.25b) - ^_:,j,fc_.)-M -"'*M ,- a J | f c- 1 (4.25c) -1 ~ » - 5 J . f c - l ^ , 4 , 4 , = K - ^ r 1 — ^ (4.25d) Ax , — 4X 1 2"'"" " 2 " ' T Z Substituting (4.25) into (4.14) eliminates Hh from the discrete linear system. 4.3.3 The Discrete Linear System This discretization can be viewed as an extension of Yee's method suit-able for discontinuous coefficients. The difference operators are centered, con-servative, and second-order accurate, as seen in [53, 54]. Throughout, a con-sistent, compact discretization of the operators grad, div, and curl is used; 84 the corresponding discrete operators can be denoted grad' 1 , div'1, curl j , and curlg. Notice that two discrete operators are needed for curl ; curl^ that op-erates on discrete vector fields defined normal to cell faces (like Ah in (4.12a-c)) and curl he acts on discrete vector fields defined to be tangential to cell edges (like Hh in (4.12g-i)). With these definitions, the identities curl Jgrad h<bh = 0h, (4.26a) d iv h cur l£H h = 0h, and (4.26b) grad hdiv h A h - curl £curl )Ah = AhAh (4.26c) follow immediately (cf. [53, 54, 64, 65, 100]). The discrete identities (4.26) are analogues of vector calculus identities. For a more thorough treatment of the finite-volume discretizations of standard vector calculus operators in (4.26), consult [64, 65]. The equations (4.10) in discrete form are cur l h e H h - (zw)- 1grad V + Jh = Jhe, (4.27a) divhjh = divhJhe, (4.27b) /j.^h - div h A h = 0h, (4.27c) jh + dh ,Ah + g r a d = Qh ( 4 2 7 d ) c u r l ) A h - iupheHh = Jhm. (4.27e) The boundary conditions (3.32) are included in the definitions of the discrete operators in (4.26) and (4.27). The operators Sh, phc, and phe are discrete approximations of the operators a and p respectively; notice two operators are needed for /z, one at cell centers, the other at cell edges. 85 After eliminating the discrete fields Hh, jh, and 4>H in (4.27), the final linear system can be formally represented as Ahxh = bh (4.28a) where A := divhdh div hdh grad hj £ j := curl^(^)- 1 curl j -grad h (^)- 1 div f c , xh := UJ (4.28b) (4.28c) (4.28d) and bh represents the appropriate source and boundary data. The operator £ ^ in (4.28c) has a 19-point stencil; the distinction between this discretization (in which Hh and iph have to be eliminated) and a direct discretization of the operator curl (yu-1 curl (•)) — grad (/x_1div (•)) is that (the value of p on cell edges) is naturally given by an arithmetic average and not a harmonic average. With p = fi0, (4.26c) implies that in (4.28d) simplifies to a scalar multiple of the discrete vector Laplacian operator Ah, so , (-itupo)-1^ - a h -dh(gvad)h \ A = \ K 1 . (4.29) \^ ( d i v ) ^ (div) ^ (grad )^ Now, applying div'1 to (4.27a) and simplifying with (4.27b) and (4.26b) gives dxvhgradhiph = 0 \ 86 Moreover, the boundary conditions (3.32c) and the discrete equation (4.27a) imply that the discrete normal derivative of tph vanishes on the boundary. Setting, say, 1^,1,1 = 0 at a single point (or applying the global constraint (3.32d)) determines that iph = 0h throughout the grid. It follows, then, that this scheme yields a discrete divergence-free vector potential, namely divhAh = Oh. (4.30) Given the exact solution (A, <p) of (3.31), the stabilization term grad V> in (4.10a) vanishes. With this scheme, (4.30) holds also, so the discrete sys-tem is indeed stabilized by this term and not penalized (e.g. [21]). This is not the case with nodal finite element methods as in [33, 69]. When an approxi-mate solution Ah fails to satisfy (4.30), the stabilization term may grow when the permeability p varies over a few orders of magnitude, or div^A^ grows in attempting to keep the discrete approximation of / u - 1 d i v A approximately constant across interfaces [69]. 4.4 Grid Effects For general physical problems, the material parameters a, p, and e are tensor fields that are functions of spatial coordinates (and possibly of frequency as well). An inherent limitation of the discretization (4.28) lies in the decision to approximate the material coefficients by functions that are piecewise constant on a tensor-product grid (cf. (4.9)). The geometry of the earth is not brick-like for realistic geophysical problems; that is, interfaces do 87 not naturally align with a simple grid. This crude model is chosen to simplify the inverse problem within which the forward-modeling problem (3.1) is embedded. For inverse problems, re-covering the material coefficient a is typically the goal, so the locations of interfaces between distinct media are not known a priori. Given that many forward solves are needed and that any physical data is sparse and contami-nated by noise, it is not reasonable to expect to be able to resolve the interfaces to arbitrary accuracy. This motivates the use of a simple model of the media like (4.9) in which the isotropic material coefficients are piecewise constant on a Cartesian tensor-product grid (4.3). More complicated finite-element discretizations could give accurate solutions for arbitrary geometries more ef-fectively but would require more a priori information about the material coef-ficients. Even highly accurate finite-volume or finite-difference discretizations are obtainable on simple Cartesian grids, but the interfaces within grid cells are assumed a priori (e.g. [74, 113]). For high-frequency electromagnetic simulations, staircase approxima-tions of interfaces are avoided (e.g. [104, Chap. 10]. In such applications, careful resolution of the geometry is crucial because errors computed will prop-agate throughout the solution [104]. This is to be expected for a hyperbolic P D E model where noise in the approximate solutions can propagate along characteristics. At low frequencies, however, Maxwell's equations are primar-ily parabolic in character, so noise in solutions tends to get smoothed out or damped by the naturally diffusive flow. Thus, the coarse approximations cho-88 sen for the material coefficients do not upset the accuracy of solutions obtained provided that the frequencies are moderate in magnitude. 4.4.1 Model Problem Analysis To what degree do solutions of Maxwell's equations change when the material coefficients are coarsely resolved as in the numerical discretizations (4.28)? To answer this question, consider first the model diffusion problem in n dimensions, -d iv (a grad cb) = f (xefl), (4.31a) with homogeneous boundary conditions 0 = 0 (xedflD), (4.31b) n • (agrad0) = 0 (x e dflN). (4.31c) In (4.31b-c), the boundary dfl = dflo U dfl^ is partitioned into two disjoint portions dflr> and dflw on which the Dirichlet and Neumann boundary con-ditions apply respectively. Inhomogeneous boundary conditions can be han-dled also (see [21, Section IV.1]) but the homogeneous boundary conditions (4.31b-c) are studied here for simplicity. Similarly, assume that the diffusion coefficient a is a positive scalar field satisfying 0 < < a(x) < cr m a x (a: G fi). The extension to the case where a is an n x n tensor field requires that a is symmetric and positive-definite with smallest eigenvalue uniformly bounded away from 0 independent of x. 89 A variational formulation of the boundary-value problem (4.31) is to find <f> G HlE(Q) such that (a g r a d 0, g r a d ip)n = [ fipdV Vip G HE(fl), (4.32a) Jn where the space of test functions is HlE{9) := {ip G ft(grad; fi) : ip = 0 on dClD} . (4.32b) Under standard assumptions on / and cr, the solution cp is known to be in HE(Q); moreover, the gradient g rad0 is also bounded (see [21, 46, 89]). Consider also a perturbed problem where the data are identical except for the coefficient cr; that is, for some h > 0, given a slightly perturbed coeffi-cient a, the perturbed boundary-value problem is -d iv (5 g r a d 0) = / ( x e f i ) , (4.33) with the associated homogeneous boundary conditions (4.31b,c) for cp. The corresponding weak problem is to find cp G HE(fl) such that (a g r a d 0, g r a d ip) = [ fipdV Vip G HE(Cl). (4.34a) ^ 'to Jn The coefficient a satisfies 0 < crm i n < a(x) < a m a x (x G fi), and cr approximates a in the sense that there is a set Sh C fi of measure 0(h) (0 < ft < 1) on which a — a is supported. That is, a = cr almost everywhere, except on a small volume near material interfaces where a fails 90 Figure 4.4: Cross-sections of the coefficient cr, the perturbed coefficient a, and the support Sh of the perturbation a — a. to resolve the interface (see Figure 4.4). With the preceding assumptions, assuming me&s(Sh) — 0(h) is the n-dimensional volume of Sh, then II" - ojl^ < K := cr m a x - <jmin and (4.35a) \\a-o-\\p<(Kmeas(Sh))r=0(hr) (1 < p < oo). (4.35b) In particular, although the discrepancy \\a — <J||p vanishes as h j 0 in standard £ p -norms , the difference is 0(1) in the £oo-norm. Unfortunately, deriving a bound on ||</» — c6||p is not straightforward in general. There is a general result bounding | | g r a d (cp — c^>) 112 -P r o p o s i t i o n 4.1. Let 4> and (p be solutions of the problems (4.32) and (4.34) respectively. Then, grad (4> - cp) Proof. With a bit of algebraic manipulation, subtracting the equations (4.32) and (4.34) yields ^ a g r a d ( 0 — 0 ) , g r a d ^ J ^ = — ((o — o) g r a d 0, g r a d i p ) n . (4.36) Substituting ip = cp — cp into the left-hand side of (4.36) gives f a\grad(cP-cP)\22dV Jn > a g r a d (cp - cp) The same substitution ip = cp — cp into the right-hand side of (4.36) yields ( ( a - a) g r a d 0 , g r a d (0-0)) < / | g r a d 0 | (a - cr)grad ( 0 - 0 ) Jn < | | g r a d 0 | | o o 11(5-a)grad ( 0 - 0 ) dV < I | g r a d 0 | | oo IKcr — cr;||2 g r a d ( 0 - 0 ) by applying Holder's inequality twice. Combining the preceding inequalities, 1 11 g r a d ( 0 < • | | g r a d 0 | | c o- - a < — | |grad0 | | O O J r r^(meas(5 ' 1 ) )3 ^min = 0(h*). • In some circumstances, the field g r a d 0 may be more relevant than 0 itself. For instance, in electrostatics, 0 would be a potential field whereas g r a d 0 would be an electric field that is measurable by physical means. The bound on | | g r a d (cp — 0)||2 in Proposition 4.1 is valuable in such cases, giving a quantitative bound on the measure of difference that staircasing makes to the computed gradients. In one dimension, it is straightforward to derive a stronger bound on | |0 - 0||oo directly. 92 Proposition 4.2. Consider the PDE problems (4.32) and (4.34) in one di-mension on the domain Cl = (0, L) for some L > 0. Assume that <f> and 0 are solutions of (4.32) with the boundary conditions 0(0) = 0'(L) = 0 (4.37) with similar boundary conditions for <f> in (4.34). Then, M-<l>U = o(h). Proof. With the boundary conditions (4.37), the exact solutions to (4.32) and (4.34) can be derived directly in one dimension, namely r F 4>(x) = / ^ Jo ° ( f l di, di, (4.38a) (4.38b) where F(i) := JLf(8)ds. Then, F(i) di < 4 - f i*(fl-crcoi-Isolde ^min 70 ll^ llc < a - o |! (4.38c) • 93 As a direct consequence of Proposition 4.2, similar bounds in the C\- and £2 _ n o r m s hold. Corollary 4.3. In one dimension, under the same assumptions as in Propo-sition 4.2, \\4> - (b\\i = 0(h) and ||0 - <f>\\2 = 0(h). Proof. The desired bounds follow directly from applying Holder's inequality on a bounded domain. In particular, The main feature that makes Proposition 4.2 possible is that the Green's functions appropriate to both the unperturbed and perturbed problems have a particular form; specifically, the solution (p in (4.38a-b) is can be written in the form where the expression ^ ( / j x ; ^ ) does not depend on a. It is reasonable to conjecture an extension of Proposition 4.2 beyond one dimension; that is, assuming the fundamental solutions of the perturbed and unperturbed P D E problems have similar structure to (4.39), the discrepancy pointwise is • (4.39) 94 and an 0(h) bound for \\(p — 4>\\oo results (as well as Corollary 4.3). However, even in one dimension, with boundary conditions other than (4.37), the term T within the Green's function depends on a (or a), i.e., T — T(o\/;x;£), making the proof slightly more complicated. 4.4.2 Numerical Study In the following tests, a reference coefficient a and a family of perturbed coefficients {crfc}£Li are prescribed in the fixed domain O — [—1, l ] 3 . For a fixed right-hand side / , associated with each of these structures are the boundary-value problems (4.31) and (4.33) (with the substitutions a = in each case). Within this section, to avoid confusion of exact and numerical solutions of perturbed and unperturbed PDEs, a discrete approximation of a continuous function 0 on a grid as described in Section 4.1 is denoted $ rather than the customary cph. Thus, <fr approximates 0, where — div (a grad <p) = / and $fc approximates <f>k, where —div (cr*, grad 4>k) = f (k = 1,... ra). Solving the preceding ra + 1 numerical problems and computing the relative discrepancies » * ; - * » ' and " g ^ * ' - * * 1 1 ' (4.40) in suitable discrete norms gives some estimation of the effect of the family of perturbations {011}™=! o n reference solution <fi. It is assumed that the dis-95 cretization is sufficiently accurate that the numerical discrepancies |<&fc — <&||p provide reasonable information about the analytic discrepancies \\<pk — (b\\p. The cell-centered finite-volume discretization coinciding with the (2,2) block of Ah in (4.28b) is applied in the experiments with a homogeneous Neu-mann boundary condition (i.e. dflo = 0 in (4.31b) and dfl^ = dfl in (4.31c)). The right-hand sides vectors fh are generated randomly. The resulting linear systems are solved with Dendy's B O X M G multigrid solver [35, 36]. For the purpose of comparison, discrete £<x>-, £\- and ^2-norms for cell-centered grid functions are respectively given by = max l$i,-fcl, (4.41a) Nx Ny Nz 1=1 j=l k=l / Nx Ny Nz \ 3 \i=i j=i k=i j Two particular reference coefficients are used in the experiments. The first is a model of a diffusive half-space o-air z> 0, anS{x,y,z) := { (4.42a) C g r o u n d Z < 0, for some positive scalars crajr and o"groun<1. For experiments with a = <THS, the coefficient in the air is fixed as crair = 10~7S/m. This test case is not difficult to solve numerically, although the jumps in a may be quite large. The simple geometry suffices to see the effect of moving the interface. 96 C H S Figure 4.5: Cross-sections of aHs and crE for perturbation experiments. For a more complicated geometry, consider also the coefficient erg de-scribing a diffusive ellipsoid buried in a half-space. In this model, oE(x,y,z) := < cra i r z > 0, ^cond S + S + ^ < 1. ( 4 4 2 b ) "ground £ <C 0, for some positive constants a C O N H, a, 6, c, and z0 with z 0 > c. The values cr a i r = 10 _ 8 S/m, (7 g r o u n d = 10 _ 3 S/m, a = b = 0.5, and c = z 0 — 0.25 are fixed in experiments. The diffusion coefficient <7E is particularly useful for a more critical evaluation of grid effects because the Cartesian grid cannot resolve the curved ellipsoid accurately. Cross-sections of the models C H S and erg are illustrated in Figure 4.5. Based on the reference coefficients (4.42), a first family of perturbations 1 to consider is based on vertical shifts of the reference model by individual grid 97 Aground N o r m $ - $1 $ _ $ 2 <"> $ _ $ 3 $ - $4 Rate II " H o c 1.451x10" 1 2.440x10" 1 2.440X10"1 2.902X10-1 0.574 10- 4 II-111 4.337x10-2 8.953x10" 2 1.346X10-1 1.788X10-1 1.033 I H I 2 3.803x10" 2 7.569x10" 2 1.126X10-1 1.492 xlO" 1 0.989 II ' I I 0 0 1.112x10-1 1.868x10" l 1.942X10-1 2.495X10"1 0.613 10- 3 II • 111 4.283x10" 2 8.340x10" 2 1.250X10"1 1.648X10-1 0.969 I I - l b 4.072x10- 2 8.207x10" 2 1.256X10-1 1.686X10"1 1.020 II • lloo 2.553x10- 1 4.091x10- 1 4.491X10"1 4.995X10-1 0.560 10- 2 II-111 5.810x10" 2 1.043x10" 1 1.491X10-1 1.920X10-1 0.855 II ' 2 6.204x10- 2 1.108x10" 1 1.555X10-1 1.983X10-1 0.837 II • lloo 1.562x10" 1 3.036x10" - l 3.535X10-1 3.820X10"1 0.782 10- 1 l l - l l i 3.319x10-2 7.126x10" -2 1.104X10-1 1.480X10-1 1.092 II • II2 3.474x10" 2 7.007x10" -2 1.051X10-1 1.379X10-1 1.005 II "lloo 2.292x10- 1 3.573x10" - l 5.050X10-1 5.210X10-1 0.651 10° II-111 8.344x10" 2 1.598x10" • 1 2.299X10"1 3.051X10-1 0.932 M I 2 7.582x10" 2 1.435x10" -1 2.055X10-1 2.709X10-1 0.916 Aground N o r m V * - v $ 1 V * - v $ 2 V $ - V $ 3 V * - V $ 4 Rate II ' l loo 6.852x10" 1 8.764x10- 1 8.765X10"1 8.765X10-1 0.252 10- 4 II'111 3.371x10" 2 6.312x10-2 9.094xl0"2 1.182X10-1 0.904 II • II2 1.187x10" 1 1.887x10" 1 2.399X10-1 2.815X10-1 0.644 II • l loo 6.057x10" l 7.630x10" 1 8.469X10"1 8.469X10-1 0.293 10- 3 3.524x10" 2 6.888x10- 2 1.008X10"1 1.321X10-1 0.959 II - 2 1.176x10" 1 1.911x10" 1 2.437X10-1 2.881X10-1 0.670 II • l loo 6.487x10" 1 8.427x10" l 8.427X10-1 8.427X10-1 0.268 10- 2 3.682x10- 2 6.905x10" 2 9.844xl0-2 1.269X10-1 0.898 M I 2 1.233x10" 1 1.954x10- 1 2.461X10"1 2.901X10-1 0.637 II ' l l o o 6.320x10- 1 9.176x10" l 9.176X10-1 9.176X10"1 0.382 10- 1 M i l 3.548x10-2 6.699x10- 2 9.637xl0-2 1.238X10"1 0.909 II ' l b 1.190x10- 1 1.903x10" 1 2.435 xlO"1 2.854X10-1 0.653 II • l loo 6.164x10- l 7.705x10- l 7.705X10-1 7.705X10-1 0.229 10° II-111 3.545x10-2 6.660x10" 2 9.627xl0"2 1.245X10-1 0.908 II • II2 1.234x10" 1 1.921x10- 1 2.452X10"1 2.889X10-1 0.626 Table 4.1: Half-space coefficient crHs'- results w: th vertical shifts 98 cells of height h. That is, ak(x,y,z) = o(x,y,z - kh) (fc = 1,...... wi). (4-43) Tables 4.1 and 4.2 present results of applying this family of perturbations (4.43) to both O"HS and erg on a grid of 72 x 72 x 72 cells. The rates of convergence are computed in each row are found by averaging fe|||||#-Ifc||| - h | # - # i | m l { k > l ) across each row (with a similar definition for the gradients). The second set of perturbed coefficients is generated by resolving the reference model on successively coarser grids. That is, given a resolved on a fine grid of spacing h, ak is obtained by resolving a on a grid (k + 1) times as coarse (i.e. of spacing (k + l)h) and injecting the values back into the original fine mesh (see Figure 4.6). The solutions $ and $ k are then computed and compared on the same fine grid of 72 x 72 x 72 cells. Table 4.3 gives results of using this family of perturbations with the ellipsoid model erg. a = <JE a = a\ a = 02 Figure 4.6: Cross-sections of ellipsoid coefficient O"E resolved on fine and coarse grids. The experiments of Tables 4.1 and 4.2 support the hypothesis that a perturbation a —a of 0(h) support induces a perturbation (p — (pof 0(h) in the 99 ^"ground Norm $ - $ 1 $ _ $ 2 $ _ $ 3 3> - $ 4 Rate II - H o c 1.620x10- l 2.986x10 - l 3.101x10-1 4.645x10-1 0.744 io- 2 3.458x10" 2 6.932x10 -2 1.106X10-1 1.561x10-1 1.049 II -112 4.507x10- 2 8.541x10 -2 1.267X10"1 1.692x10-1 0.939 II • lloo 2.241x10" 1 3.171x10 - l 4-OlOxlO- 1 4.536x10-1 0.513 10- 1 II • 111 3.120x10" 2 6.644x10 -2 1.050X10-1 1.420x10-1 1.096 II " 112 4.416x10- 2 8.414x10 -2 1.239X10-1 1.630x10-1 0.937 II • lloo 1.630x10- 1 2.688x10 -1 3.323x10-1 3.673x10-1 0.652 10° II • 11 x 3.673x10" 2 7.358x10 -2 1.097x10-1 1.455x10-1 0.997 II • l b 4.546x10" 2 8.837x10 -2 1.286x10-1 1.663x10-1 0.947 II • lloo 1.446x10" 1 2.498x10 - l 3.145x10-1 3.617x10-1 0.719 101 II-111 6.289x10- 2 1.312x10 -1 2.008x10-1 2.706x10-1 1.057 II • II2 5.667x10- 2 1.168x10 -1 1.780x10-1 2.392x10-1 1.041 II • lloo 1.548x10" l 2.905x10 - l 3.420x10-1 3.823x10-1 0.760 102 3.868x10- 2 7.754x10" -2 1.159x10-1 1.544x10-1 1.000 II • II2 4.444x10" 2 8.672x10 -2 1.276x10-1 1.681x10-1 0.962 Aground Norm vs- v $ 1 2 V $ - V $ 3 V $ - V $ 4 Rate II • lloo 6.481x10" 1 8.785x10" i 8.785x10-1 8.785x10-1 0.312 10- 2 II-111 3.799x10-2 7.093x10- 2 l . O l l x l O - i 1.307x10-1 0.894 II • l b 1.278x10-1 1.992x10- 1 2.498x10-1 2.939x10"! 0.617 II • lloo 6.982x10- 1 8.153x10- 1 8.153x10-1 8.636x10-1 0.173 10- 1 3.678x10" 2 6.962x10" 2 1.001x10-1 1.303x10-1 0.915 M I 2 1.243x10- 1 1.988x10" 1 2.508x10"! 2.972x10-1 0.648 II "l loo 6.181x10" l 7.459x10" l 8.274x10-1 8.274x10-1 0.249 10° M l i 3.890x10" 2 7.391x10" 2 1.044X10-1 1.322x10-1 0.903 • II2 1.257x10- 1 2.023x10- 1 2.549X10-1 2.964X10-1 0.650 II • lloo 6.288x10- l 8.107x10" 1 8.357X10-1 8.513x10-1 0.281 101 M l i 3.624x10-2 6.980x10- 2 1.003X10-1 1.296x10-1 0.931 II ' 2 1.200x10-1 1.973x10" 1 2.489x10-1 2.948x10-1 0.676 II ' lloo 6.278x10- l 9.572x10" l 9.572x10-1 9.572x10-! 0.432 102 IJ • jl! 4.064x10" 2 7.712x10" 2 1.103x10-1 1.424x10-1 0.913 II • » 2 1.258x10- 1 2.044x10- 1 2.582x10-1 3.040x10-1 0.664 Table 4.2: Ellipsoid coefficient aE- results with vertical shifts 100 Ci- and /Vnorms. Further, the / V n o r m field grad (<p — 0) seems to decrease with h at a rate close to O(h^) and clearly less than 0(h), consistent with the prediction of Proposition 4.1. There is a marked difference in the convergence of the £oo-norms of <& — v e r s u s the other norms; this is indicative of the fact that the difference is very strongly concentrated near the set Sh and rapidly decays away from there. Measurement of convergence rates in Table 4.3 is not as straightforward where the ellipsoidal model aE is resolved on a family of successively coarser grids. However, it is reassuring to notice that the relative differences in the numerical solutions are extremely small even with the coarse resolution of the ellipsoid's curvature. The results shown here demonstrate that staircasing is not as problematic in diffusive problems as it is in dispersive ones; such is the case for Maxwell's equations considered at sufficiently low frequencies also. 101 Aground N o r m $ - $1 $ _ $ 2 $ _ $3 Rate II • H o c 1.331x10" ti 2.148x10" ti 2.190x10- ti 0.572 i o - 2 II-111 3.321x10- 6 6.001x10" 6 4.060x10- 6 0.518 II -112 1.779x10" 6 3.221x10" 6 2.240x10- 6 0.533 II ' IIco 1.400x10- ti 1.775x10" ti 1.860x10- ti 0.300 i o - 1 II-Hi 1.705x10" 7 3.447x10" 7 4.141x10- 7 0.911 I I ' 2 2.107x10" 7 3.567x10" 7 4.603x10- 7 0.735 II • H o c 2.232x10- ti 4.279x10" ti 4.287x10- ti 0.294 10° M l i 3.410x10" 7 7.221x10" 7 1.078x10" 6 1.074 II • l b 4.398x10" 7 7.566x10" 7 1.073x10" 6 0.726 II ' l l o o 1.999x10" ti 2.476x10" ti 2.715x10" ti 0.326 10 1 II-Hi 2.270x10" 7 4.901x10" 7 7.111x10" 7 1.260 II • II2 3.252x10" 7 5.337x10--7 7.311x10" 7 0.850 II • l l o o 2.147x10" ti 2.506x10-- t i 3.443x10" ti 0.767 10 2 II-111 2.473x10-7 5.977x10--7 9.732x10" 7 1.065 II • II2 3.632x10" 7 6.437x10--7 9.497x10" •7 0.797 Aground N o r m 2 3 4 Rate II ' l l o o 4.985x10- ti 6.917x10" ti 7.087x10- ti 0.396 10- 2 II • 111 3.906x10-7 6.271x10" 7 8.363x10- 7 0.688 II • II2 1.104x10-6 1.495x10- 6 1.858x10" 6 0.456 II • l l o o 5.803x10" ti 6.446x10- ti 6.721x10- ti 0.143 10- 1 l l - l l i 4.631x10-7 7.136x10- 7 9.450x10- 7 0.637 1.300x10- 6 1.730x10- 6 2.101x10- 6 0.425 II • l l o o 6.574x10- ti 1.107x10" b 1.107x10" i> 0.613 10° ' 1 4.770x10- 7 7.398x10- 7 9.938x10" 7 0.651 1.306x10" 6 1.784x10" 6 2.178x10- 6 0.458 II ' l l o o 5.781x10" ti 6.386x10" ti 6.655x10" ti 0.136 ro1 M l i 4.456x10- 7 7.039x10" 7 9.403x10" 7 0.670 II ' II2 1.300x10" 6 1.737x10- 6 2.122x10" 6 0.432 II • l l o o 6.546x10" ti 6.566x10- ti 7.256x10" ti 0.049 102 M l i 4.550x10" 7 7.127x10- 7 9.808x10" 7 0.673 II < II2 1.309x10" 6 1.733x10- 6 2.153x10- 6 0.429 Table 4.3: Ellipsoid coefficient <JB: results with fine & coarse resolutions 102 Chapter 5 Construction and Analysis of Iterative Solvers Since the frequency-domain problem (3.1) needs to be solved many times for various frequencies and sources within a single geophysical inverse problem, the major objective of this work is to construct fast solvers for the forward-modeling problem. A major obstacle in achieving this goal is the ill-conditioning of the corresponding discretizations due to the near-vanishing conductivity cr in a large portion of the domain. Another obstacle is the strong coupling between solution components as a result of the double-curl P D E operator in (3.1a) Finally, the presence of strong discontinuities in the coefficients a (and possibly p) contributes also to ill-conditioning in the linear systems that can inhibit solver performance. The P D E formulation and dis-cretization in preceding chapters deal with these concerns in part. It remains to find robust algorithms for the solution of the linear algebraic equations. Among iterative methods for solving linear systems of algebraic equa-tions, two particular families stand out for this application: preconditioned 103 Krylov-subspace methods [10, 49, 93] and multigrid methods [19, 18, 22, 56]. In considering these strategies, a particular multigrid preconditioner approx-imating the dominant diagonal blocks of the coefficient matrix is a natural choice for an iterative solver. After development of a Fourier analysis in the spirit of [17, 24, 38] to determine bounds of the spectrum of the preconditioned system under certain assumptions, some numerical experiments demonstrate the efficacy of this solver. 5.1 Iterative Methods With the analytic reformulation of (3.1) as the boundary-value problem (3.31), the discretization (4.28) derived in Chapter 4 yields a set of equations Ax = b (5.1) describing the discrete fields (Ah, <ph) over a three-dimensional domain. For a typical problem using a grid of 40 x 40 x 40 cells, there are over 250,000 complex unknowns to determine. In the absence of unusual structure in the sparse coefficient matrix, problems of this size rule out direct solvers, so iterative solvers are required for (5.1). 5.1.1 Preconditioned Krylov-Subspace Methods With a lexicographic ordering to the unknowns (Ah, <ph), the basic block 104 structure of the coefficient matrix A in (5.1) is \ (5.2) An — iLuBn —IU1B12 A j B21 A22 J In the general case, An is a matrix representation of the real discrete difference operator £ j = c u r l £ ( / 4 ) - 1 c u r l J - grad' l (^)- 1 div' 1 as in (4.28c). The operator has a 19-point stencil that simplifies to a 7-point stencil for the vector Laplacian operator in the case that p = / J 0 , as in (4.29). The block A22 represents the difference operator d iv h ah grad h as in (4.28), while the matrix blocks Bn, B12, and B21 are discrete representations of lower-order operators. Given that the diagonal blocks An and A22 correspond to discretizations of second-order differential operators, for a sufficiently fine grid, these blocks are dominant. This dominance is strongly dependent on the frequency u>; for high frequencies, the grid spacing h has to be quite small to ensure dominance of An in the top block row of (5.2). With suitable boundary conditions and scaling, the coefficient matrix A in (5.2) is complex symmetric, but non-Hermitian and indefinite. As such, the well-known method of conjugate gradients and the associated theory cannot be applied to solve the system of linear equations [49, 61, 93]. Even so, non-Hermitian Krylov-subspace iterations such as BiCGStab or G M R E S are known to perform well for a wide range of applications with suitable precondition-ers (i.e. matrices M such that M approximates A - 1 ) to improve convergence behavior (e.g. [38, 87, 62]). 105 As a simple example of a preconditioner for (5.1), consider a single step of SSOR (Symmetric Successive Over-Relaxation [34, 10]). Like other precon-ditioners based on a single step of a classical, stationary iterative method, an SSOR preconditioner is simple to implement, requiring only knowledge of the splitting A = D - L - U, (5.3a) where L and U are strictly lower- and upper-triangular matrices and D is a diagonal matrix. With the splitting (5.3a), the SSOR preconditioner is Ms(a) := (D-aU)-1(aL + (l-a)D)(D-aL)-1(aU + (l-a)D), (5.3b) where a G (0,2) is the relaxation parameter1 [10]. For certain second-order elliptic P D E problems, a choice of the relaxation parameter aopt exists that reduces the spectral condition number of cond(Ms(o:optM) optimally; however, practical computation of a o p t requires information about the spectrum of A that is prohibitively expensive [48]. If the forward-modeling problem (5.1) must be solved repeatedly for many different frequencies u>, the preconditioner Ms(a) must be computed for each UJ since both D and U in (5.3a) depend on UJ. The block structure of A suggests a different preconditioner. Given the diagonal dominance of the blocks A n and A22 in (5.2), the matrix M, E \ A22 J (5.4) :The relaxation parameter is denoted with a here rather than the customary symbol UJ to avoid confusion with the frequency. 106 is an approximation of A'1 to leading order. The frequency dependent blocks IUJBII, iuB22, and VJOB\2 are ignored in defining ME in (5.4). If numerous systems of the form (5.1) need to be solved for different frequencies, since ME is independent of UJ, then only a single L£/-factorization of ME need be computed and stored. Of course, for large three-dimensional problems, computing ME exactly by LU-factorization of An is almost as much work as doing computing an L U -factorization of A itself. As such, some approximation of ME is sought. A well-known method for doing this is through applying ILU (Incomplete LU) factorizations to the diagonal blocks of A (e.g. [93, Chap. 10]). That is, rather than finding the exact LU-factorizations of An and A22, approximate factors £11^11 — -An and L22U22 - A 22 are computed. The triangular factors Lu and Uu (i = 1,2) are computed as in the usual LU-factorization, and either fill elements below a prescribed drop tolerance are rejected (ILU with a threshold) or all fill elements are disregarded (ILU(O)). The block-ILU preconditioner Mr := y L22U22 j (5.5) turns out to be an effective approximation of ME that in turn approximates A^1 for sufficiently low frequencies. By ignoring the off-diagonal blocks of A, computing M j requires much less work and storage than a similar ILU-factorization of A. 107 The efficacy of the ILU factorization approximating ME comes in part from some additional structure; in particular, both of the blocks of ME come from the discretization of second-order elliptic differential equations. For such problems, multigrid methods have been extensively studied and devel-oped that exploit the smoothing properties of the underlying elliptic operators [19, 18, 56, 113]. As an alternative to applying a preconditioned Krylov-subspace iteration to solve the linear system (5.1), a full multigrid method could be developed to solve the system (4.28). That is, a multigrid method could approximate the inverse of A directly in (5.1) rather than choosing a pre-conditioner M and applying a Krylov-subspace iteration. Such a full multigrid method could be derived through traditional geometric arguments from the discrete operator in (4.28), or algebraically directly from the matrix represen-tation A in (5.1). The geometric-multigrid approach is somewhat difficult for systems of PDEs discretized on a staggered grid; development of the inter-grid transfer operators is fairly straightforward, but distributive smoothers in three dimensions require more effort [18, 17]. The algebraic-multigrid approach is more amenable to a black-box implementation and is adaptable to unstruc-tured grids as well; however, in algebraic-multigrid solvers, the construction of a hierarchy of coarse grids is not obvious and the theoretical details are often more complicated than for geometric multigrid [92]. The results of an experiment applying an A M G solver are documented in Section 5.4. A sensible middle ground adopted here is to apply a multigrid method as a preconditioner for a Krylov-subspace iteration [114, 9, 32, 62, 108]. In 108 this particular practical multigrid preconditioner MM approximates the action of ME in (5.4) by a multigrid method. That is, MM := Mu \ V M-22 (5.6) where Mu represents the action of a single multigrid V(2,l)-cycle approximat-ing the action of A^1 (i=l,2). For a homogeneous permeability /i = po, the block An in (5.2) reduces to a 3 x 3 block-diagonal matrix. In that case, the block preconditioner ME is a 4 x 4 block-diagonal matrix, as are Mi in (5.5) and MM in (5.6). The multigrid preconditioner, then, is MM •--M, 22 M, 33 V (5.7) M 4 4 / The diagonal blocks of MM in (5.7) represent the action of multigrid methods for solving four uncoupled P D E problems in the scalar fields Ax, Ay, Az, and <p (specifically, three Poisson problems and one inhomogeneous diffusion problem). Constructing MM in (5.7) is simpler than in (5.6) because MN requires building a multigrid solver for a 3 x 3 coupled P D E system discretized on a staggered grid. The preconditioner MM in (5.7) is analyzed in the next section. For the practical implementation in the experiments of Section 5.4, a black-box multigrid solver—in particular, the code B O X M G of Dendy [35, 36]—is applied. 109 5.2 Analytical Framework To better understand the preconditioner ME, consider the system (3.31) and its discretization (5.1) under simplifying assumptions that ensure that the eigenvectors of all the blocks of A are Fourier modes. That is, for discrete wavenumbers 0 < a, (3,7 < ft — 1, the grid function V>Q/37 given by ^ : = e(2^aifc).e(2^jfc).e(2^rfch). ij,k = 0,l,...n-l (5.8) should be an eigenvector of all the blocks of A in (5.1). This is true under the assumptions: •1. The boundary conditions are periodic. 2. The grid is uniform. 3. The material properties cr, /i, and e are all constant. Under the first assumption, upon discretization, the indexing is done in arithmetic modulo ft. For instance, grid points on, say, the top boundary are identified with opposite points on the bottom boundary, and centers of cells touching the top boundary are considered as neighbors of centers of cells touching the bottom boundary. The scalar grid functions AXH, AYH, AZH, and <ph each consist of N := ft3 unknowns, the matrix blocks are square of dimension N x N, and the resulting system to be solved has dimension 4iV x AN. The first assumption seems unreasonable in light of the boundary condi-tions (3.32) actually implemented. However, for a class of realistic geophysical 110 problems, the fields axe expected to vanish or decay in the absence of source currents, and the sources are frequently of compact support. Thus, for a suf-ficiently large domain O, the effect of boundary data on the fields inside the domain is often negligible. In such cases, the assumption of periodicity in the following analysis is particularly reasonable and does not compromise the model. Further, as in Section 2.2.4, make the quasistatic assumption a = a and make a slight regularizing assumption a <— max(<7, <Tair), (5.9) with < 7 A I R > 0 small (on the order of 10 - 9 to 10 - 6 S/m). The regularizing assumption (5.9) ensures a ^ 0 everywhere (which would not be the case ordinarily in the air). These assumptions simplify the presentation without loss of essence in the obtained theoretical results. In the experiments reported in Section 5.4, these simplifying assumptions on a are relaxed as are most other assumptions of the present section. 5.2.1 Discrete Operators on a Periodic Grid ' With the preceding assumptions, the matrix (4.28) reduces to (4.29) since the permeability LL is constant. Furthermore, since the conductivity a is also constant, the block matrix ah reduces to al. The entire coefficient matrix 111 (5.10) A can be written as ((iup0)-lAh + aIN S aD^ (iu)ii0)-lAh + aIN uDhy (iufi0)-1Ah + aIN aDhz -<j(Dhxy -o-{D*y -o-{Dhzy oAh) where the matrices D%, Dy, and Dz correspond to discretizations of first-order derivative operators with periodic boundary conditions in (5.10). That is, assuming standard lexicographic ordering of the unknowns describing the grid functions, let V D-1 h -1 1 G C (5.11) v - 1 i / be the matrix for the periodic primitive backwards difference operator in one dimension. Then, Dt := /„ In ® D~ G C NxN Dhy := In ® D~ ® In G C Dt N * N , and D - ® l n ® l n e C N x N , where ® denotes the conventional Kronecker product (e.g. [34, p. 274]). Like-wise, the matrices — (£)£)*, — (Dy)*, and —(Dt)* in the fourth block row of (5.10) correspond to forward difference operators. Finally, the discretization Ah of the Laplacian operator A with periodic boundary conditions satisfies Ah = -Dhx{Dhxy - Dhy{Dhyy - Dhz(Dhzy. 112 The primary advantage of considering the discretization (4.28) on a periodic lattice is that the discrete difference operators in (5.10) share common eigenvectors. Explicitly, for a, ft, 7 — 0, . . . , n — 1, the discrete grid functions tpa^y defined by : = e ( 2 T O i h ) , e ( ^ h ) . e ( ih r 7 * fc ) . (ij,k = 0,1,...n-1) are common eigenvectors of the difference operators D%, Dy, D^, Ah, and of the identity matrix IN- The corresponding eigenvalues are dxaf* := h~l (1 - e - { 2 a 7 T h ) l ) for£>£, (5.12a) dya^ := h'1 (1 - e - ( 2 f 3 n h ) l ) forDj , (5.12b) d°fr := hr1 (1 - e - { 2 ^ h ) l ) for (5.12c) := __L (Sin2(a7r/i) + sin 2(/?7r/i) + sin2(77r/i)) for A \ (5.12d) In particular, the nonzero eigenvalues 5 Q / ? 7 of A' 1 satisfy 0 < ^ < - ^ < l ( « , / ? , 7 = 0 , . . . , n - l ; a 2 + / ? 2 + 7 2 > 0 ) (5.13) for n > 2. The upper bound of 1/16 on —l/8al31 is obtained for a very coarse mesh when n = 2; as n increases, this upper bound decreases. With the eigenvalues (5.12) of the blocks of A in (5.10), determine now the eigenvalues and singular values of the preconditioned iteration matrix M~lA, where M := Ah Ah \ oAh (5.14) 113 The matrix M is real and, up to a scalar multiple, consists of the dominant diagonal blocks of the matrix A, excluding the Helmholtz shift al^ in the first three block rows. Note that the solution u of Au — b with periodic boundary conditions is determined only up to a constant. Thus, M in (5.14) has a zero eigenvalue of multiplicity four (corresponding to four distinct constant grid functions for Axh, Ayh, Azh, and <bh). However, as in [38], if So denotes the space spanned by the constant grid function ip000 € and SQ denotes the orthogonal com-plement of So in C ^ , then A ^ x is invertible. Further, the operator ( A f e | 5 x ) _ 1 extends to an operator (A*1)'1 defined over all of by defining it to be 0 on So- Thus, it makes sense to consider the matrix M~lA restricted to a subspace S that excludes the 4-dimensional constant null-space of M , as in (5.13). Denote this matrix by M~lA\s. 5.3 Spectral Bounds The goal now is to derive bounds on the eigenvalue range and the I2-condition number of M~lA\s; this is possible due to the fact that the blocks of A in (5.10) all share common eigenvectors. The eigenvalues and singular values of this matrix are determined by considering the eigenvalues of a 4 x 4 matrix with the same structure, as the following lemma indicates. L e m m a 5.1. Let A G C(MN>X<MN> be a block matrix with blocks Atj G Cnxn (i,j = 1,... ,m). Let the blocks A\j all share a common eigenvector v G C " with corresponding eigenvalues G C (z,j = l,...,m). Furthermore let 114 B = (ccij) E C m x m have the eigenvalues ctij as its elements, and let A be an eigenvalue of B with an eigenvector u E C m . Then A is also an eigenvalue of A with an eigenvector u®v. Proof. Although this lemma is standard, the proof is included for complete-ness. By definition, AijV = otijV (i, j = 1, 2, . . . , m) Bu — Xu. Thus, the iih block of A(u <g> v) is m m 3=1 3=1 = XuiV = A [ i t (8) w ] i . • Next, find the eigenvalues and singular values of the 4 x 4 matrix cor-responding to M - 1 A | s in Fourier space. Introduce the parameter X:=a;/xcr>0 (5.15) that has dimensions of L~2, where L is a unit of length. L e m m a 5.2. Let dx, dy, dz EC be prescribed scalars such that 5 :— —dxdx* — dydy* — dzdz* < 0. Given o, x > 0, define the complex 4 x 4 matrix M_1A, 115 where M := \ , and A := ad) Then, the eigenvalues X of M lA are A 2 = A 3 := 1 + 7«, o 4^ := 1 + ^ (x ~ VxU ~ 4^)) . Further, the singular values r\ of M lA are m = c 1 l + - + - V ^ + 4 ) i + 82 c 1 1 + 2 - 2 ^ + 4) <5 + zx zxrfx 5 + zx «x^y 5 + «x ^X^-2 ^—ado;* —crdy* —adz* a5 J (5.16a) (5.16b) (5.16c) (5.17a) (5.17b) (5.17c) where c = c ( - J , x ) := ( l + X2 ( l - J)) > 0. Proof. Although the derivation of (5.16) and (5.17) requires tedious algebraic manipulations, the formulae are verifiable directly by inspection. • Notice that the eigenvalues and singular values of M~lA in Lemma 5.2 depend only on the parameters 5 > 0 and x- Next, derive bounds for |A| and rj in terms of 5 and x-116 Lemma 5.3. Assume 0 < — 1/5 < JQ. Then, the eigenvalues A of M lA given in (5.16) all satisfy 0 < A m i n < A — A m a x for some positive constants that depend on x but n°t on 5. Similarly, the singular values rj given in (5.17) all satisfy max for some positive constants that depend on x but not on 5. Further, as x —> oo, A r A m i r = 0(X2) and = 0(X2) Vmin Proof. Denote C = \Jx2 + 1652 > X- From (5.16), it is apparent that \/2x , Ai — A 2 A 4 = X 452 A? (x + C + V^Vc + x) - ^ ^ x +1 l + — <52 (x + C - V^xVC + x) + ^ p x / C ^ x +1 It follows from the analytic expressions in (5.18) that IAi| > |A 2 | = |A 3 | > |A 4 | for any valid choices of x a n d 5. Examining |Ai| and |A 4 | in (5.18a,c) as 5 varies over the interval (—oo, —16], it follows that , (5.18a) (5.18b) . (5.18c) \nax ; — max {|Ai| : 5e ( - o o , - 1 6 ] J = |Ai Amin : = min | | A 4 | : 5 G (—oo, —16]} = | A 4 5=-16 <5=—16 , and (5.19a) (5.19b) 117 Thus, using the analytic expressions in (5.19), as y_ ~^ °°> + 6 + 0(^) = 0(X2). A max Ai A 4 Y2 1 ^ - + ( -256 x Amir Examining the singular values (5.17) in the same manner yields Vi>m = r]3> VA-Thus, ?7max := max {r/i : 5 G (-oo, -16]} = r7i|5=_16 , and (5.20a) ?7 m i n := min{?74 : 8 G (-oo, -16]} = ry 4 | 5 = _ 1 6 . (5.20b) Hence, the bounds in (5.20) yield • Lemma 5.3 indicates that the ratio of the largest to smallest eigenvalues significantly underestimates the ratio of the largest to smallest singular values (which is the ^-condition number) in the case that x 1S large. In Figure 5.1, the curves plotted are based on the above estimates for the eigenvalue range and the ^-condition number of M~lA as a function of x- Also, Figure 5.1 shows the actual values (in Fourier space) for a given instance n = 104. With the bounds established for the 4 x 4 case, Lemma 5.1 gives the eigenvalues and singular values of the full coefficient matrix M~lA\s-118 — I?ma*/>?mi„ ° (n = 104) "max min x A , n a I t / A m l n (n = 10") / o / , 'o ,' 10 10 Figure 5.1: Eigenvalue range A m a x / A m i n and the condition number rjmax/r]mm as a function of x, based on the formulae in Lemma 5.4; also plotted are actual values ( A m a x / A m i n and r}m^/nmin, respectively) for n = 104. L e m m a 5.4. The eigenvalues of M 1A\s are given by A f 7 = 1 + X + Vx(x ~ 4<5^z)) , X A f 7 = A f : = l + ^ i , A 4 1 • (x ~ VX{X - 4 < J ° ^ ) ) (5.21a) (5.21b) (5.21c) (a, (3,7 = 0, . . . , n — 1; a2 + L32 + 7 2 > 0). The singular values of M 1A\s are given by a/37 1 + — - + - x / c a / 3 7 ( c Q ^ + 4) V2 = % • = 1 + r,01^ — 74 C<*P1 \ 1 + — - v / c a ^ ( c Q ^ + 4) (5.22a) (5.22b) (5.22c) where ca^ = c ( - ^ , X ) := (-^ W) ( l + xM 1 -^) ) > °-Proof. Apply Lemma 5.1 using the scalar eigenvalues and singular values of 119 the 4 x 4 matrix from Lemma 5.2 along with the eigenvalues (5.12) of the difference operators constituting (5.10). • Combining the results of the preceding lemmata gives the final result. Theorem 5.5. The eigenvalues and singular values of M~lA\s are bounded above and below independent of the grid size h. Further, the eigenvalue range Amax/Amin a n ^ ^ e condition number r}max/r]min are bounded independent of the grid, but grow quadratically with the parameter x '•= co/ia; i.e. ^ = 0(X2) and = 0(x2). ^min Vmin The practical implication of Theorem 5.5 is that the number of iter-ations required for convergence (up to a set positive tolerance) of a Krylov method for (5.10) preconditioned by the real parts of the leading blocks is expected to be small and independent of the grid size. As the frequency u increases in a range where copo >> 1, the number of iterations is expected to grow like u?. The full inversions of the Laplacians on the main block diagonal of M may be prohibitively expensive as a preconditioner. Instead, consider apply-ing just one multigrid cycle. Denote by B the approximate inverse for M\s obtained in this way. It is well known (e.g. [56, 113]) that for a W-cycle us-ing standard grid transfer and coarse grid operators and a sufficiently good smoothing rate (e.g. two red-black Gauss-Seidel sweeps before and after the coarse grid correction), there is a constant c independent of the grid size such 120 that \\I - BM\S\\ < c < 1. This implies that the condition number of BM\s is bounded. Writing cond(BA) < c o n d ( £ M | s ) cond(M\slA) and applying Theorem 5.5 now yields similar results for the one multigrid cycle preconditioner. Corollary 5.6. The condition number of BA\s, the problem matrix precondi-tioned by one multigrid cycle, is bounded independent of the grid size h. This condition number increases quadratically with the parameter x '•= ujpo~. As a result of the preceding corollary, it is hoped that at low- to mid-range frequencies UJ, using a few multigrid W- or F-cycles to approximate the inver-sion of the diagonal blocks yields a robust preconditioner for the linear system (5.1). 5.4 Numerical Experiments In the experiments that follow, the parameters e = eo = 8.8542 x 10 - 1 2 F / m and the permittivity p = po = An x 10 - 7 H / m are fixed constants. The grids used are either uniform or exponentially-widening (as described in [6]). The domain fl is the cube [—1, l ] 3 and the homogeneous boundary conditions (3.32) are imposed on dQ, as described in Chapter 3. Although use of the boundary conditions (3.32) requires a domain fl that is sufficiently 121 large, the experiments here focus on the performance of iterative solvers based on multigrid algorithms; this is unrelated to the specific boundary conditions implemented [18]. In comparing the performances of BiCGStab, QMR, and G M R E S [10] in solving the linear systems of [6], BiCGStab performs significantly better; hence, consistent with the experiments of [6, 53, 54], preconditioned BiCGStab is chosen as the Krylov-subspace solver for the linear systems in the exper-iments that follow. As a convergence criterion for all BiCGStab iterations, a reduction of the relative residual to within a tolerance of 10~7 is required. Results from all experiments are given in terms of iteration counts. Several preconditioning strategies are compared for performance. 1. Exact solves of diagonal blocks ( M E ) : The real parts of the diagonal blocks, namely diag (Ah + u2el, Ah + oj2eI, Ah + u2el, (div) V ( g r a d )h) , are inverted exactly in each iteration by iterating a multigrid solver to convergence. This differs slightly from (5.4) due to the nonvanishing Helmholtz shift u>2e. For the range of frequencies considered, the dif-ference between this preconditioning strategy and the one analyzed in Section 5.3 is negligible. 2. Multigrid preconditioning (MM)' Single V(2,l)-cycles are applied for each of the diagonal blocks to approximate the effect of M E as in (5.7). For the (1,1), (2,2), and (3,3) blocks of MM, point red-black Gauss-Seidel 122 smoothing is used if the grid is uniform. Otherwise, alternating-plane smoothing is used. For the (4,4) block, alternating-plane smoothing is used in all cases. 3. ILU preconditioning (Mr): An incomplete-LU factorization is computed using the real part of the matrix M E as in (5.5). A drop tolerance 10 - 2 is applied. 4. SSOR preconditioning (Ms): The value a = 1.0 used for the over-relaxation parameter in (5.3b). The multigrid solver used is Dendy's code B O X M G (a black-box multi-grid solver, see [35, 36]). For Poisson or simple Helmholtz problems on uniform grids, point red-black Gauss-Seidel smoothing cheaply inverts the (1,1), (2, 2), and (3, 3) blocks. For the inhomogeneous diffusion operator in the (4,4) block and for problems on non-uniform grids, alternating-plane smoothing is required to achieve reasonable multigrid performance, as observed in [1, 35, 36]. This gives a slight cause for concern due to the greater cost of plane relaxations. However, for large, highly inhomogeneous problems, the multigrid precondi-tioner applied with alternating-plane smoothing yields a robust solver. In terms of computing costs (floating point operations or wall-clock time), the SSOR preconditioner Ms is cheapest. The ILU preconditioner Mr is slightly more expensive per iteration (depending on the level of fill) and requires some cost to set up, but is still fairly cheap. When on uniform grids and point red-black Gauss-Seidel smoothing is used on the Laplacian blocks, the multigrid preconditioner MM takes roughly four times as long per iteration 123 as the preconditioner M j ; for nonuniform grids, alternating-plane relaxation is needed, so each iteration of MM is roughly 8.5 times as long as one of Mr, ^cond Aground Figure 5.2: Cross-section of CTB ( a conducting block in a half-space). As a first model, define the conductivity O"B in fl by " a i r , Z > 0] "cond, N < | , ! b l | < | , - | < ^ < 0 ; ( 5- 2 3) " g r o u n d , otherwise. This models a conducting block of conductivity <Jcond embedded in a half-space. For the second test, consider the conductivity O"E describing an ellipsoid buried underground (cf. (4.42b)). This constitutes a more difficult problem because the material interfaces do not align with the grid. The conductivity a g r o u n d of the earth is c r g r 0 u n d = 10~3 S/m and the conductivity < 7 a i r of the air is " a i r = 10"' S/m (which regularizes a as in (5.9)). Thus, the parameters c r c o n d and UJ are varied to compare performances of various preconditioners. The applied electromagnetic source term Je in (3.31) is either an electric or a magnetic dipole source as in [53, 110], while Jm = 0. Consider first varying UJ over the moderate range [10°, 104] Hz, and 124 BiCGStab iterations 0"cond #of ijj = 10 uHz OJ = 10 2 Hz OJ = 10 4 Hz (S/m) cells ME MM Mi M s ME MM Mi Ms ME MM Mi Ms 10 3 1 3 7 19 3 3 7 19 6 6 11 20 20 3 2 2 13 34 2 2 13 34 5 5 22 35 io- 2 30 3 2 2 20 44 3 3 20 44 5 5 34 50 40 3 3 3 27 58 3 3 27 61 5 5 46 70 50 3 3 3 33 72 3 3 33 79 5 5 61 86 10 a 3 3 9 25 3 3 9 25 6 6 12 25 20 3 3 3 18 44 3 3 19 44 5 5 35 48 10° 30 3 3 3 30 66 3 3 30 50 5 5 56 63 40 3 3 3 41 73 3 3 41 73 5 5 66 82 50 3 3 3 44 100 3 3 44 100 5 5 100 111 10 3 3 3 9 20 3 3 9 25 6 6 12 25 20 3 3 3 20 44 3 3 20 44 6 6 37 48 10 2 30 3 3 3 30 67 3 3 27 59 6 6 55 67 40 3 3 3 40 73 3 3 40 73 6 6 76 86 50 3 3 3 49 95 3 3 48 96 6 6 97 109 Table 5.1: a = <7B and electric dipole source on uniform grids. Each iteration of MM is roughly 4 times as long as a single iteration of Mj. varying crccmd over the range [10 - 2,102] S/m. These ranges translate to the parameter y_ m (5.15) lying in the interval [47T X 10 _ 9,47T X 10_ 1] m - 2 . The jump discontinuities in cr for this range of crcond are between 5 and 9 orders of magnitude. For an electric dipole source, the results when discretized on uniform and nonuniform grids are tabulated in Tables 5.1 and 5.2 respectively. Tables 5.3 and 5.4 present results from similar trials on uniform and nonuni-form grids respectively with a divergence-free magnetic dipole source term instead. Table 5.5 presents results of a similar set of trials with the model CE discretized on a set of uniform grids with an electric dipole source term. All of the trials summarized in Tables 5.1-5.5 clearly show evidence of grid-independent rates of convergence of the Krylov method when precondi-tioned by ME or MM, as opposed to MT or Ms. In spite of the seemingly strong 125 BiCGStab iterations " c o n d # o f w = 10 u Hz 10 2 Hz w = 10 4 Hz (S/m) cells M E M M Mi Ms ME M / M S ME MM Mi M S 30 3 1 3 18 47 3 :s 18 49 5 5 34 52 io-2 40 3 4 4 27 59 4 4 26 64 7 7 40 68 50 3 2 2 28 76 2 2 28 71 4 4 48 84 30 3 3 3 26 65 3 3 26 62 5 5 55 76 10° 40 3 4 4 31 76 4 4 34 78 7 7 76 91 50 3 2 2 43 101 2 2 43 102 4 4 98 139 30 3 3 3 25 65 3 3 25 64 6 6 53 76 10 2 40 3 4 4 31 79 4 4 34 80 7 7 76 102 50 3 2 2 43 105 2 2 43 102 5 5 100 132 Table 5.2: a = <TB and electric dipole source on nonuniform grids. Each iteration of MM is roughly 8.5 times as long as a single iteration of M j . restrictions of the Fourier analysis of Section 5.2, the convergence behavior of the exact block inversion ME measures up to theoretical predictions. That is, the analysis succeeds in capturing the essential local behavior of discrete op-erators for physically reasonable problems unhindered by the additional com-plications of nonuniform grids, of different boundary conditions, and, most significantly, of highly discontinuous conductivity models. The results also show virtually identical performance in terms of it-eration counts for the exact preconditioner ME versus the single multigrid V(2,l)-cycle MM- This observation suggests that solving the diagonal blocks exactly is unnecessary, at least within this frequency range. In particular, single V(2,l)-cycles work as effectively as the preconditioner ME, with consid-erable savings in computational work. Another interesting aspect of these results is that the iteration counts for all the preconditioners do not appear to depend strongly on u or acon(x-Looking at Figure 5.1, the spread of eigenvalues and the ^-condition number 126 BiCGStab iterations ^cond # o f U) = 10 u Hz u> = 10 2 Hz a> = 10 4 Hz (S/m) cells ME Mj Ms M E MM Mi Ms ME MM Mi Ms 10 3 1 4 6 13 4 4 6 13 4 4 6 13 20 3 4 4 13 25 4 4 13 25 4 4 15 25 io-2 30 3 4 4 19 32 4 4 19 32 4 4 21 32 40 3 4 4 27 42 4 4 27 42 4 4 32 41 50 3 4 4 32 51 4 4 32 53 4 4 36 53 10 3 4 4 6 13 4 4 6 13 4 4 6 13 20 3 4 4 13 25 4 4 13 25 4 4 15 25 10° 30 3 4 4 19 32 4 4 19 32 4 4 21 32 40 3 4 4 27 43 4 4 27 41 4 4 32 43 50 3 4 4 32 51 4 4 32 51 4 4 37 48 10 3 5 5 7 13 5 5 7 13 5 5 7 13 20 3 5 5 13 26 5 5 13 26 5 5 15 26 10 2 30 3 4 4 19 32 4 4 19 32 4 4 22 32 40 3 5 5 27 42 5 5 27 42 5 5 32 43 50 3 5 5 31 52 5 5 31 50 5 5 37 56 Table 5.3: a = crB and magnetic dipole source on uniform grids. Each iteration of MM is roughly 4 times as long as a single iteration of MT. BiCGStab iterations ^cond # o f U) 10 u Hz w = 10 2 Hz UJ = 10 4 Hz (S/m) cells ME Mi Ms ME MM Mi Ms M E MM Mi Ms 30 3 1 3 23 41 3 3 23 42 3 3 25 42 io-2 40 3 3 3 30 45 3 3 30 50 3 3 36 51 50 3 2 2 37 61 2 2 37 60 2 2 47 62 30 3 3 3 23 42 3 3 23 41 3 3 26 41 10° 40 3 3 3 30 47 3 3 30 44 3 3 37 48 50 3 3 3 37 60 3 3 37 65 3 3 47 65 30 3 3 3 21 41 3 3 21 41 4 4 27 48 10 2 40 3 4 4 30 54 4 4 30 49 4 4 37 49 50 3 3 3 37 57 3 3 37 57 4 4 47 59 Table 5.4: a = <7B and magnetic dipole source on nonuniform grids. Each iteration of MM is roughly 8.5 times as long as a single iteration of Mr. 127 BiCGStab iterations *^cond #of UJ = 10uHz UJ = lO^Hz UJ = 104Hz (S/m) cells MB MM M 7 ME MM ME M M Mi 303 2 4 20 2 4 20 4 8 72 i o - 2 403 2 4 30 2 4 30 4 6 107 503 2 4 34 2 4 34 4 6 146 303 2 4 18 4 4 18 9 6 64 10° 403 2 4 25 2 4 24 5 6 98 503 2 4 29 3 4 33 5 5 133 303 3 4 19 4 5 23 9 10 59 102 403 2 4 24 2 4 29 8 9 92 503 2 4 33 2 4 33 9 8 133 Table 5.5: a = CE arid electric dipole source on uniform grids. Each iteration of MM is roughly 4 times as long as a single iteration of M / . is quite flat until x — 101 m - 2 . Once x gets larger, however, the quadratic dependence of both quantities on x becomes apparent and we can see the condition number getting far worse. Since x is at most 47r x 10 _ 1 m - 2 in the trials so far, the iterations do not change much as crconcj and UJ are varied within this range. BiCGStab iterations #of UJ (Hz) cells 101 10^ 103 104 10b 10° 303 3 3 3 6 12 98 MM 403 3 3 3 6 13 116 503 3 3 3 6 12 128 303 30 27 31 55 166 642 Mj 403 40 40 42 76 210 1180 503 46 48 51 97 273 1551 Table 5.6: a — erg over a range of frequencies. To better observe the effect of increasing the parameter x-, fix the un-derlying conductivity model with cr c o n d = 102 S/m and vary the frequency UJ over a larger range using an electric dipole source. The multigrid pre-128 conditioner MM and the ILU preconditioner Mr are compared for problems discretized on uniform grids and the results are summarized in Table 5.6. Ta-ble 5.6 demonstrates the degradation of both preconditioning strategies MM and Mi with increasing frequency. As in the other trials, iteration counts for the ILU-preconditioned system grow as the grid is refined, as opposed to the grid-independent iteration counts observed for the multigrid-preconditioned system. At the highest frequency to = 106 Hz, x ~ 47r m - 2 . Moreover, the performance of both preconditioners worsens significantly in increasing to from 105 Hz to 106 Hz. This is consistent with the sharp rise observed in r/max/^ rnin in Figure 5.1 near x — 10 m - 2 . Thus, the observations for higher frequencies agree well with the predictions of the Fourier analysis. This example is set over a short length scale, but the same quantitative results apply in terms of iteration counts on larger domains. The CPU-time required for each iteration2 of the multigrid precondi-tioner is roughly 4 times more than that for the ILU preconditioner in case of a uniform grid, and roughly 8.5 times more in the case of the exponentially widening grid (due to the additional cost of alternating-plane relaxations). Thus, using the ILU preconditioner Mr actually results in a faster algorithm than using the multigrid preconditioner MM for coarse discretizations. But for fine discretizations (i.e., all cases recorded in Table 5.6) the multigrid pre-conditioner is significantly more efficient. Obviously, there is a crossover point beyond which the multigrid preconditioner is more efficient, since it achieves 2 The iteration times are compared by running trials with the preconditioners MM and Mi on various machines. Precise computational work estimates for B O X M G are in [35, 36]. 129 convergence to a fixed tolerance within a constant number of iterations inde-pendent of the grid spacing h. In this example, the crossover point occurs for grid sizes that are well within practical range. A more specific implementa-tion of the multigrid cycle could improve the efficiency of the preconditioner further, and lead to an earlier crossover point. As a final test, Figure 5.3 presents convergence curves comparing the performance of MM-preconditioned BiCGStab with an algebraic multigrid code3. The test problem uses the ellipsoidal model coefficient <TE from (4.42b) discretized on a grid of 30 x 30 x 30 cells. The code B O X M G is applied as a black-box solver to the dominant diagonal blocks of A in (5.1) while the other code applies an algebraic multigrid method to the full system. The re-sults for MM are plotted for several frequencies, while those for the algebraic solver (denoted by A M G in the Figure 5.3) are only plotted for one frequency; with the algebraic multigrid method used, the convergence curves are virtually identical for the distinct frequencies considered. Although restricted to a small test case due to implementation details, the results summarized in Figure 5.3 are important in demonstrating the favorable convergence behavior for MM in spite of the cruder approximation inherent in this preconditioning strategy. The multigrid-preconditioning strategy turns out to be competitive with the more intensive full A M G solver in spite of the added simplicity (indeed, for this comparison, the wall-clock time for A M G is nearly ten times greater). This emphasizes the utility in the approach adopted throughout this research of deriving preconditioners from the structure of the P D E model as opposed 3The AMG code was generously provided by Yair Shapira. 130 to more brute force solvers. Convergence of A M G and MM 1 i ! \ ~ ~ ^ - ^ ' ' V. \ 1 1 1 1 1 A M G (UJ = 103) 1 0.01 MM (UJ = 101) —h— CM L i \ ' \ \ MM (UJ = 102) • X -\ 1 » o 0.0001 i\ '* ^ \ M M (UJ = 103) * bX) O \ x. X le-06 \ x X X . . 1 ¥ | , , 0 2 i 6 8 10 12 14 16 18 Iterations Figure 5.3: Comparison of A M G solver and MM-131 Chapter 6 Conclusions 6.1 Research Summary Starting from the mathematical foundations of electromagnetic theory, this research addresses various concerns relevant to low-frequency electromag-netic simulation in geophysical or medical inverse problems. Drawing on the experience of practitioners of computational fluid dynamics, a P D E descrip-tion suitable for moderate-frequency regimes is derived. With an appropriate continuous formulation, a finite-volume discretization is presented that pos-sesses the conservation properties of the established Yee scheme and provides a sensible framework for respecting subtle interface conditions. A study of the consequences of staircase approximations to complex geometries in diffusive problems highlights the feasibility of Cartesian grids widely applied in prac-tice. Finally, a family of preconditioners for the iterative solution of the result-ing systems of equations is constructed, analyzed, implemented, and tested, yielding reliable solvers for large-scale three-dimensional problems. 132 An underlying contention of this work is that the construction of reli-able numerical schemes and algorithms for the solution of PDEs is inextricably linked to the selection of physically meaningful and mathematically consistent forms of the associated continuous problems. The assumptions underlying the P D E model are critically examined from this perspective in Chapter 3. For the ranges of parameters considered, second-order boundary-value PDEs in the magnetic field H are ruled out due to the strong discontinuities in the con-ductivity cr. To circumvent difficulties caused by the nontrivial kernel of the dominant curl operator, the electric field is split with a Helmholtz decomposi-tion into divergence-free and curl-free components expressed in terms of vector and scalar potential variables A and (b. The resulting indefinite P D E system is subsequently stabilized by subtracting a vanishing correction grad (p~l div A). Applying techniques established in the study of incompressible fluid flow, an extra differentiation yields a diffusion P D E for the scalar potential <p resem-bling a pressure-Poisson equation. The resulting diffusion equation replaces the gauge condition, yielding a P D E operator in Cartesian coordinates with much weaker coupling between solution components than the original P D E in the standard variables E or H. The new P D E problem is completed with two extra boundary conditions: one is needed to specify the Helmholtz decomposi-tion and the other is required due to the additional differentiation. The logical progression outlined in Chapter 3 provides the analytical foundation on which the subsequent discretization and iterative methods stand. With the P D E (3.31) expressed in terms of scalar and vector poten-133 tials, the finite-volume discretization (4.28) is derived as in [53, 54]. As with the Yee scheme, the important conservation properties of the discretization (4.28) result from careful interlacing of field components on a grid staggered in space. In addition, the finite-volume arguments used to derive the discrete equations (4.28) provide a logical framework for the discrete representations of discontinuous material coefficients. As a result, discrete solutions approximate necessary interface conditions (3.34) at the junctions of distinct physical me-dia. Furthermore, the coefficient matrix of the linear system stemming from the discretization (4.28) inherits block-diagonal dominance from the ellipticity and structure of the P D E operator in (3.31) expressed in Cartesian coordi-nates. The structure of this linear system is amenable to the preconditioners described in Chapter 5. Discretization schemes widely used in practice for inverse problems are based on the assumption that the isotropic, inhomogeneous material coeffi-cients are piecewise constant on grid cells. With standard Cartesian tensor-product grids, this supposition limits the resolution of interfaces between phys-ical media to staircase approximations aligned with the faces of the cells of the grid. For straightforward scalar inhomogeneous diffusion problems, Proposi-tion 4.1 states that if an interface between media is perturbed by 0(h), then the gradient of the solution of a perturbed diffusion problem differs from the gradient of a solution of an unperturbed diffusion problem by O(h^) in the £2-riorm. Furthermore, in one dimension, the solution of the diffusion problem itself is perturbed by 0(h) in the £\-, C2-, and ^ co-norms; this is particularly 134 surprising given that a perturbation of 0(h) of the interface generally results in a perturbation of 0(1) in the £oo-norm of diffusion coefficient. While ex-tension of the one-dimensional result to higher dimensions is not obvious, the numerical study of Section 4.4.2 suggests the same result holds in three di-mensions. Thus, the implicit assumption of insensitivity of P D E solutions to staircase approximations of material interfaces is shown to be justifiable within the range of diffusive problems considered. With the numerical scheme (4.28) and the corresponding sparse sys-tem of linear algebraic equations established, various preconditioned Krylov-subspace methods are considered for the iterative solution of the equations. A preconditioning strategy that stands out is one based on isolating the domi-nant diagonal blocks of the coefficient matrix. These blocks correspond to dis-cretizations of uncoupled, scalar, elliptic PDEs with well-established solution techniques. Two robust preconditioners emerge that approximate inversions of the diagonal blocks, one based on ILU-factorization and another using a multigrid method. The reliability of these block-diagonal preconditioners fol-lows directly from the properties of the P D E operator at moderate frequencies. The Fourier analysis of Section 5.2 supports the preceding assertion. Under the assumptions outlined in Section 5.2, the coefficient matrix of the discrete linear system has a block structure wherein all blocks share Fourier modes as eigenvectors. As a result, Theorem 5.5 provides explicit bounds for the eigenvalue range and the ^-condition number independent of the grid spac-ing. Moreover, these bounds are shown to grow quadratically with x — wfM>~, 135 so degradation in the performance of block diagonal preconditioning strategies is expected as the frequency grows. In spite of the extreme limitations of the Fourier analysis, experiments in Section 5.4 strongly support the theoretical predictions for a variety of cases. Thus, for sources having compact support, the action near the boundary is negligible, so the local analysis captures the essential properties of the discrete operator well even in the presence of strong spatial inhomogeneities. Further, the block diagonal solver is observed to be competitive with an algebraic multigrid solver for the full discrete system, em-phasizing the strength of the simpler block-diagonal preconditioning strategy. 6.2 Future Research The formulation of the boundary-value problem (3.31) in terms of the vector and scalar potentials (A, cp) is presented in Chapter 3 with the explicit assumption of existence of solutions. For Maxwell's equations expressed in the usual field variables E and H, establishing the existence and regularity of solutions in distinct physical contexts remains an area of active research in the mathematical community (e.g. [4, 97]). Complications within the geophysical context include the rough nature of the inhomogeneous coefficients and solu-tions and the selection of appropriate boundary conditions. For the purpose of applied computations within this work, proceeding with experimentation based on the assumption of well-posedness is justifiable. However, for the sake of mathematical consistency, the existence of solutions in appropriate spaces should be established. 136 Another mathematical question to address is the convergence of the finite-volume scheme (4.28). The derivation provided in Chapter 4 and the results of [54, 53] show that the method is consistent and second-order accu-rate; however, as is well-known in the numerical analysis community, given consistency of a numerical scheme, stability (in whatever appropriate sense) is generally required to conclude convergence. It is possible that searching for a relationship with established techniques using finite elements will point the way to a convergence proof (e.g., as in [81]). The analysis and experiments of Section 4.4 raise a few interesting ques-tions. It remains to extend Proposition 4.2 beyond the simple O D E case for one dimension; that is, is it true that a perturbation of 0(h) in the resolu-tion of the interface between distinct media results in a perturbation of 0(h) in the solution of the inhomogeneous diffusion problem in the £oo-norm? A similar theoretical result for more general elliptic P D E problems (including time-harmonic Maxwell's equations) would also be of practical utility. For more general boundary conditions and geometries in higher dimensions, more sophisticated analytical tools are needed to push the result up. Approximate boundary conditions in computational electromagnetics constitute a broad area of research. A n assumption implicit in many of the experiments and the analysis of this thesis is that the accuracy of the chosen boundary conditions for Maxwell's equations is not crucial for the family of applications studied herein. Although outside the scope of this work, it might be worthwhile to examine the suitability of this assumption with a numeri-137 cal study comparing the solutions of Maxwell's equations with P M L , P M C , absorbing, or other boundary conditions. There are some examples of such comparisons at high frequencies [31, 104, 105] but not at lower frequencies of interest here. Another possible extension of the current work would be in the devel-opment of a time-domain code for Maxwell's equations. Some geophysical surveys rely on time-domain models of macroscopic electromagnetic behavior rather than time-harmonic frequency-domain models. For instance, in tran-sient electromagnetic surveys, a D C source current in a steady state is abruptly shut off. The practical reason for shutting off the source in this manner is to allow receivers to measure the secondary electromagnetic fields due to cur-rents induced in conductors without measuring the primary field caused by the source (that may be several orders of magnitude larger [82]). For transient electromagnetic phenomena and problems with discontinuous coefficients, the frequency-domain P D E models, spatial discretizations (see also [6, 53]) and solvers studied herein could be adapted for the time domain. For parabolic problems of interest, it is likely that implicit time-stepping is needed, unlike the explicit time-stepping used with the F D - T D method [104, 100, 112]. That being the case, the solvers described in Chapter 5 are of great utility; see [11, 70, 63] among other efforts. 138 Bibliography [1] R.E . Alcouffe, A. Brandt, J .E. Dendy, and J.W. Painter. The multi-grid methods for the diffusion equation with strongly discontinuous coeffi-cients. SIAM J. Sci. Stat. Comput, 2:430-454, 1981. [2] D. Alumbaugh and F. Morrison. Electromagnetic conductivity imaging with an iterative born inversion. IEEE Trans, on geoscience and remote sensing, 31:758-762, 1993. [3] D. Alumbaugh, G. Newman, L. Prevost, and J.N. Shadid. Three di-mensional wideband electromagnetic modeling on a massively parallel computers. Radio Science, 31:1-23, 1996. [4] H. Ammari, A. Buffa, and J.-C. Nedelec. A justification of eddy currents model for the Maxwell's equations. SIAM J. Appl. Math., 60:1805-1823, 2000. [5] D. Aruliah and U. Ascher. Multigrid preconditioning for time-harmonic Maxwell's equations, submitted, 2000. [6] D. Aruliah, U. Ascher, E . Haber, and D. Oldenburg. A method for the forward modelling of 3D electromagnetic quasi-static problems. M3AS, 11:1-21, 2001. [7] U. Ascher, R. Mattheij, and R. Russell. Numerical Solution of Boundary Value Problems for Ordinary Differential Equations. SIAM, 1995. [8] U. Ascher and L. Petzold. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations. SIAM, Philadelphia, PA, 1998. [9] O. Axelsson and P. S. Vassilevski. A survey of multilevel preconditioned iterative methods. BIT, 29:769-793, 1989. 139 [10] R. Barrett, M . Berry, T . F . Chan, J. Demmel, , J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, and H. Van der Vorst. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. SIAM, 1994. [11] R. Beck, P. Deuflhard, R. Hiptmair, R. Hoppe, and B. Wohlmuth. Adap-tive multilevel methods for edge element discretizations of Maxwell's equations. Surveys Math. Indust., 8:271-312, 1999. [12] J.-P. Berenger. A perfectly matched layer for the absorption of electro-magnetic waves. J. Comput. Phys., 114:185-200, 1994. [13] O. Biro and K. Preis. On the use of the magnetic vector potential in the finite element analysis of three-dimensional eddy currents. IEEE Trans., MAG, 25:3145-3159, 1989. [14] M . Born. Principles of Optics. Pergamon Press, 1975. [15] A. Bossavit. Solving Maxwell's equations in a closed cavity, and the question of "spurious modes". IEEE Trans., MAG, 26:702-705, 1990. [16] A. Bossavit. Computational electromagnetism. Variational formulation, complementarity, edge elements. Academic Press, 1998. [17] A. Brandt. Guide to multigrid development. In W. Hackbusch and U. Trottenberg, editors, Multigrid Methods, pages 220-312. Springer-Verlag, 1982. [18] A. Brandt. Multigrid techniques: 1984 Guide with applications to fluid Dynamics, 1984. [19] A. Brandt and N . Dinar. Multigrid solution to elliptic flow problems. In Numerical Methods for PDE's, pages 53-147. Academic Press, 1979. [20] F. Brezzi, J. Douglas Jr., and L.D. Marini. Two families of mixed finite elements for second order elliptic problems. Numer. Math., 47:217-235, 1985. [21] F. Brezzi and M . Fortin. Mixed and Hybrid Finite Element Methods. Springer-Verlag, 1991. 140 [22] William L. Briggs. A Multigrid Tutorial. SIAM, 1987. [23] G. Carey and T. Oden. Finite Elements, a second course. Prentice-Hall, 1983. [24] T . F . Chan and H.C. Elman. Fourier analysis of iterative methods for elliptic boundary value problems. SIAM Review, 31:20-49, 1989. [25] M . Clemens, R. Schuhmann, U. van Rienen, and T. Weiland. Modern Krylov subspace methods in electromagnetic field computation using the finite integration theory. ACES J., 11:70-84, 1996. [26] M . Clemens and T. Weiland. Comparison of Krylov-type methods for complex linear systems applied to high-voltage problems. IEEE Trans., MAG, 34:3335-3338, 1998. [27] M . Clemens and T. Weiland. Numerical algorithms for the FDITD and F D F D simulation of slowly varying electromagnetic fields. Int. J. Num. Modelling, 12:3-22, 1999. [28] M . Clemens and T. Weiland. Transient Eddy-Current Calculation with the FI-Method. IEEE Trans., MAG, 35:1163-1166, 1999. [29] M . Costabel. A Coercive Bilinear Form for Maxwell's Equations. J. Math. Anal. Appl, 2:528-541, 1991. [30] R. Courant and D. Hilbert. Partial Differential Equations. Interscience Publishers, 1962. [31] F. El Dabaghi, K. Morgan, A .K. Parrott, and J. Periaux, editors. Ap-proximations and Numerical Methods for the Solution of Maxwell's Equa-tions. Oxford University Press, 1998. [32] W. Dahmen and A. Kunoth. Multilevel preconditioning. Numer. Math., 63:315-344, 1992. [33] N. A. Demerdash and R. Wang. Theoretical and numerical difficulties in 3D vector potential methods in finite element magnetostatic compu-tations. IEEE Trans., MAG, MAG-26T656-1658, 1990. [34] J. Demmel. Applied Numerical Linear Algebra. SIAM, 1997. 141 [35] J .E . Dendy, Jr. Black Box Multigrid. J. Comput. Phys., 48:366-386, 1982. [36] J .E . Dendy, Jr. Two multigrid methods for three-dimensional problems. with discontinuous and anisotropic coefficients. SIAM J. Sci. Stat. Comp., 8(2):673-685, Sept. 1987. [37] A. Dey and H.F. Morrison. Resistivity modeling for arbitrarily shaped three dimensional structures. Geophysics, 44:753-780, 1979. [38] H. Elman. Preconditioning for the steady-state Navier-Stokes equations with low viscosity. SIAM J. Scient. Comput, 20:1299-1316, 1999. [39] H. Elman and D. Silvester. Fast Nonsymmetric Iterations and Pre-conditioning for Navier-Stokes Equations. SIAM J. Scient. Comput, 17:33-46, 1996. [40] H.W. Engl, M . Hanke, and A. Neubauer. Regularization of Inverse Prob-lems. Kluwer, 1996. [41] B. Engquist and A. Madja. Absorbing boundary conditions for the nu-merical simulation of waves. Math. Comp., 31:1-24, 1977. [42] C. Farquharson and D. Oldenburg. Non-linear inversion using general measures of data misfit and model structure. Geophysics J., 134:213-227, 1998. [43] R.P. Feynman, R.B. Leighton, and M . Sands. The Feynman Lectures on Physics. Addison-Wesley, 1977. [44] C .A.J . Fletcher. Computational Techniques for Fluid Dynamics, vol-ume II. Springer-Verlag, 1988. [45] R. Freund and N. Nachtigal. QMR: A quasi-minimum residual method for non-Hermitian linear systems. Numer. Math., 60:315-339, 1991. [46] V. Girault and P.A. Raviart. Finite Element Methods for Navier-Stokes Equations. Springer, 1986. [47] D. Givoli. Numerical Methods for Problems in Infinite Domains. Else-vier, 1992. 142 Gene H. Golub and Charles F. Van Loan. Matrix Computations. Johns Hopkins Univ. Press, 1996. Anne Greenbaum. Iterative methods for solving linear systems. SIAM, 1997. P.M. Gresho and R.L. Sani. On pressure boundary conditions for the incompressible Navier-Stokes equations. Internat. J. Numer. Methods Fluids, 7:1111-1145, 1987. P.M. Gresho and R.L. Sani. Incompressible Flow and the Finite Element Method: Advection-Diffusion and Isothermal Laminar Flow. John Wiley & Sons Inc., 1998. E. Haber. Numerical Strategies for the Solution of Inverse Problems. PhD thesis, University of British Columbia, 1997. E. Haber and U. Ascher. Fast finite volume simulation of 3D electro-magnetic problems with highly discontinuous coefficients. SIAM J. Sci-ent. Comput., 2001. To appear. E . Haber, U. Ascher, D. Aruliah, and D. Oldenburg. Fast simulation of 3D electromagnetic using potentials. J. Comput. Phys., 163:150-171, 2000. E . Haber, U. Ascher, and D. Oldenburg. On optimization techniques for solving nonlinear inverse problems. Inverse problems, 2000. In press. W. Hackbusch. Multi-Grid methods and Applications. Springer-Verlag, 1985. T. Hagstrom, B. Alpert, L. Greengard, and S. Hariharan. Accurate Boundary Treatments for Maxwell's Equations and their Computational Complexity. 14th Annual Review of Progress in Applied Computational Electromagnetics, to appear. [58] T. Hagstrom and S. Hariharan. A formulation of asymptotic and exact boundary conditions using local operators. Appl. Num. Math., 27:403-416, 1998. 143 [59] R.F. Harrington. Field Computation by Moment Methods. MacMillan, 1968. [60] M . Hestenes. Conjugate Direction Methods in Optimization. Springer-Verlag, 1980. [61] M . Hestenes and E . Stiefel. Methods of conjugate gradients for solving linear systems. Cons. Reas. Nat. Bur. Stand., 49:409-436, 1952. [62] R. Hiptmair. Multilevel preconditioning for mixed problems in three di-mensions. PhD thesis, University of Augsburg, 1996. [63] R. Hiptmair. Multigrid method for Maxwell's equations. SIAM J. Nu-mer. Anal, 36:204-225, 1998. [64] J. M . Hyman and M . Shashkov. Natural discretizations for the diver-gence, gradient and curl on logically rectangular grids. Int. J. Computers & Math, with Appl, 33:81-104, 1997. [65] J .M. Hyman and M . Shashkov. The adjoint operators for natural dis-cretizations for the divergence, gradient and curl on logically rectangular grids. IMACS J. Applied Numerical Math., 25:413-442, 1997. [66] J .M. Hyman and M . Shashkov. Mimetic discretizations for Maxwell's equations. J. Comp. Phys., 151:881-909, 1999. [67] R.S. Ingarden and A. Jamiolkowski. Classical Electrodynamics. Elsevier, 1985. [68] J.D. Jackson. Classical Electrodynamics. John Wiley k. Sons Inc., 1975. [69] J. Jin. The Finite Element in Electromagnetics. John Wiley & Sons Inc., 1993. [70] P. Ciarlet Jr and J. Zou. Fully discrete finite element approaches for time-dependent Maxwell's equations. Numer. Math., 82:193-219, 1999. [71] A . A . Kaufman and G.V. Keller. The magnetotelluric sounding method. Elsevier, 1981. 144 [72] D.J. LaBrecque. A scalar-vector potential solution for 3D E M finite-difference modeling. In M . Oristaglio and B. Spies, editors, Interna-tional Symposium on Three-Dimensional Electromagnetics, pages 143-149. Schlumberger-Doll Research, 1995. Rolf Leis. Initial boundary value problems in mathematical physics. John Wiley & Sons Inc., 1986. R.J. LeVeque. Numerical Methods for Conservation Laws. Birkhauser, 1990. Y. Li and D.W. Oldenburg. Inversion of 3-D DC resistivity data us-ing an approximate inverse mapping. Geophysical journal international, 116:557-569, 1994. P. Lorrain, D.P. Corson, and F. Lorrain. Electromagnetic Fields and Waves. W.H. Freeman and Company, 1988. R.L. Mackie, T.R. Madden, and P.E. Wannamaker. Three dimensional magnetotelluric modeling using difference equations - theory and com-parison to integral equation solutions. Geophysics, 58:215-226, 1993. J.T. Marti. Introduction to Spbolev Spaces and the Finite Element So-lution of Elliptic Boundary Value Problems. Academic Press, 1986. C. Mattiussi. An analysis of finite volume, finite element, and finite difference methods using some concepts from algebraic topology. Journal of Computational Physics, 133:289-309, 1997. P. Monk and E. Suli. A convergence analysis of Yee's scheme on nonuni-form grids. SIAM J. Numer. Anal, 31:393-412, 1994. J.D. Moulton. Numerical implementation of nodal methods. PhD thesis, Institute of Applied Mathematics, University of British Columbia, 1996. M.N. Nabighian and J.C. Macnae. Time domain electromagnetic prospecting methods. In Electromagnetic methods in applied geophysics: applications, pages 427-520. Soc. Expl. Geophys., 1991. 145 [83] G. Newman and D. Alumbaugh. Three-dimensional massively parallel electromagnetic inversion-I. Theory. Geophysical journal international, 128:345-354, 1997. [84] G.A. Newman and D.L. Alumbaugh. Frequency-domain modelling of air-borne electromagnetic responses using staggered finite differences. Geo-phys. Prospecting, 43:1021-1042, 1995. [85] D. W. Oldenburg. One-dimensional inversion of natural source magne-totelluric measurements. Geophysics, 44:1218-1244, 1979. [86] R. L. Parker. Geophysical Inverse Theory. Princeton University Press, Princeton NJ, 1994. [87] I. Perugia and V. Simoncini. An optimal indefinite precondi-tioner for a mixed finite element method. Technical Report 1098, Publ. I.A.N. C.N.R., 1998. [88] I. Perugia, V. Simoncini, and M . Arioli. Linear algebra methods in a mixed approximation of magnetostatic problems. SIAM J. Scient. Comput, 21 (3): 1085-1101, 1999. [89] M . Renardy and R. Rogers. An Introduction to Partial Diferential Equa-tions. Springer-Verlag, 1993. [90] S.M. Richardson and A.R.H. Cornish. Solution of three-dimensional incompressible flow problems. J. Fluid Mech., 82:309-319, 1977. [91] R.D. Richtmyer and K.W. Morton. Difference Methods for Initial-Value Problems. John Wiley and Sons, Inc., 1967. [92] J. W. Ruge and K. Stiiben. Algebraic multigrid (AMG). In S. F. Mc-Cormick, editor, Multigrid Methods, pages 73-130. SIAM, 1987. [93] Y. Saad. Iterative Methods for Sparse Linear Systems. PWS Publishing Company, 1996. [94] Y. Saad and M.H. Schultz. GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comp., 7:856-869, 1986. [95] J .E . Santos and D. Sheen. On the existence and uniqueness of solu-tions to Maxwell's equations in bounded domains with applications to magnetotellurics. M3MAS, 10:615-628, 2000. [96] S. Selberherr. Analysis and Simulation of Semiconductor Devices. Springer-Verlag, 1984. [97] D. Sheen. Approximation of electromagnetic fields: Part I. Continuous problems. SIAM J. Appl. Math., 57:1716-1736, 1997. [98] J. Shen. Computational Electromagnetics using Boundary Elements. Computational Mechanics Publications, 1994. [99] D. Sidilkover and U. Ascher. A multigrid solver for the steady state Navier-Stokes equations using the pressure-Poisson formulation. Comp. Appl. Math. (SBMAC), 14:21-35, 1995. [100] T. J. Smith. Conservative modeling of 3D electromagnetic fields; part I: properties and error analysis. Geophysics, 61:1308-1318, 1996. [101] D.C. Stinson. Intermediate Mathematics of Electrmagnetics. Prentice-Hall, 1976. [102] G. Strang and G.J. Fix. An Analysis of the Finite Element Method. Prentice-Hall, 1973. [103] J.C. Strikwerda. Finite Difference Schemes and Partial Differential Equations. Wadsworth & Brooks/Cole, 1989. [104] A. Taflove. Computational Electrodynamics: the Finite-Difference Time-Domain Method. Artech House Publishers, 1995. [105] Allen Taflove, editor. Advances in computational electromagnetics: the finite-difference time-domain method. John Wiley & Sons Inc., 1998. [106] S.V. Tsynkov. Numerical Solution of Problems on Unbounded Do-mains. A review. Appl. Numer. Math., 27:465-532, 1998. [107] C R . Vogel. Numerical solution of a non-linear ill-posed problem arising in inverse scattering. Inverse Problems, 1, 1985. 147 [108] C. Vuik, J.J.I.M. van Kan, and P. Wesseling. A black box multi-grid preconditioner for second order elliptic partial differential equa-tions. In G. Bugeda E. Onate and B. Suarez, editors, Proceedings Euro-pean Congress on Computational Methods in Applied Sciences and En-gineering, ECCOMAS 2000, Barcelona, Spain, September 11-14, 2000. CIMNE, 2000. [109] T. Wang and G.W. Hohmann. A finite-difference time-domain for three-dimensional electromagnetic modelling. Geophysics, 58:797-809, 1993. [110] S.H. Ward and G.W. Hohmann. Electromagnetic theory for geophysical applications. In Electromagnetic Methods in Applied Geophysics, pages 131-311. Soc. Expl. Geophys., 1988. [Ill] C. Weaver. Mathematical Methods for Geo-Electromagnetic Induction. John Wiley & Sons Inc, Philadelphia, 1994. [112] J .T. Weaver and A . K . Agarwal. Recent developments in three dimen-sional finite difference modelling of the magnetic field in geo-electric induction. International Symposium on Three Dimensional Electromag-netics, pages 131-142, 1995. [113] P. Wesseling. An Introduction to Multigrid Methods. John Wiley & Sons Inc., 1992. [114] J. Xu and J. Qin. Some remarks on a multigrid preconditioner. SIAM J. Sci. Comput, 15:172-184, 1994. [115] K.S. Yee. Numerical solution of initial boundary value problems involv-ing Maxwell's equations in isotropic media. IEEE Trans, on antennas and propagation, 14:302-307, 1966. 148
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Fast solvers for time-harmonic Maxwell’s equations...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Fast solvers for time-harmonic Maxwell’s equations in 3D Aruliah, Dhavide Arjunan 2001
pdf
Page Metadata
Item Metadata
Title | Fast solvers for time-harmonic Maxwell’s equations in 3D |
Creator |
Aruliah, Dhavide Arjunan |
Date Issued | 2001 |
Description | The speed of iterative solvers for discretizations of partial differential equations (PDEs) is a significant bottleneck in the performance of codes designed to solve large-scale electromagnetic inverse problems. A single data inversion requires solving Maxwell's equations dozens if not hundreds of times. An inherent difficulty in geophysical contexts is that the conductivity and permeability coefficients may exhibit discontinuities spanning several orders of magnitude. Furthermore, in the air, the conductivity effectively vanishes. In standard formulations of Maxwell's equations, the curl operator that dominates the PDE operator leads to strong mixing of field components and illconditioning of linear systems resulting from standard discretizations. The primary objective of this research is to build fast iterative solvers for the forward-modeling problem associated with electromagnetic inverse problems in the frequency domain. Toward this goal, a Helmholtz decomposition of the electric field using a Coulomb gauge condition recasts the PDE problem in terms of scalar and vector potentials. The resulting indefinite system is then stabilized by addition of a vanishing term that lies in the kernel of the dominant curl operator. Finally, an extra differentiation recasts the PDE system in a diagonally-dominant form reminiscent of a "pressure-Poisson" formulation for incompressible fluid flow. The continuous PDE problem obtained is equivalent to the original Maxwell's system but has a structure that is amenable to reliable solution techniques. Using a finite-volume scheme, the PDE is discretized on a staggered grid in three dimensions. The discretization obtained possesses conservation properties typical of finite-volume methods. Furthermore, interface conditions imposed by discontinuities in the material coefficients are sensibly accounted for in deriving the discretization. Although the simple representation of the media on a Cartesian tensor-product grid uses staircase approximations of surfaces of discontinuity of the material coefficients, some analysis and a numerical study demonstrate the suitability of such coarse approximations for diffusive problems. The discretization yields a non-Hermitian sparse linear system of algebraic equations; various preconditioners for Krylov-subspace methods are described, analyzed, implemented, and tested. Of particular interest is a multigrid preconditioner that exploits both the structure of the PDE problem and the availability of well-established solvers for elliptic PDE problems (in particular, Dendy's BOXMG solver). The end result is a robust solver for the forward-modeling equations that can be incorporated within a competitive inverse problem code. |
Extent | 8185516 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-10-01 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0051685 |
URI | http://hdl.handle.net/2429/13522 |
Degree |
Doctor of Philosophy - PhD |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-714276.pdf [ 7.81MB ]
- Metadata
- JSON: 831-1.0051685.json
- JSON-LD: 831-1.0051685-ld.json
- RDF/XML (Pretty): 831-1.0051685-rdf.xml
- RDF/JSON: 831-1.0051685-rdf.json
- Turtle: 831-1.0051685-turtle.txt
- N-Triples: 831-1.0051685-rdf-ntriples.txt
- Original Record: 831-1.0051685-source.json
- Full Text
- 831-1.0051685-fulltext.txt
- Citation
- 831-1.0051685.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0051685/manifest