Quasi-minimal Residual Research Articles

A general forward solver for 3D CSEMs with multitype sources and operating environments

To determine the electromagnetic (EM) fields of different three-dimensional (3D) controlled-source electromagnetic methods (CSEMs) using the same parameters of the forward solution, by explicitly considering the commonalities, we present a general 3D forward modeling solver for CSEMs with multitype sources and operating environments. The commonality of the solver is reflected in two aspects. First, the solver is based on a frequency-domain (FD) vector Helmholtz equation for determining the scattered electric field. The different types of sources are imposed on the right-hand term of the equation, expressed as background Green’s function. Second, sources of any CSEM can be composed of electric dipole (ED) or magnetic dipole (MD) superposition. Thus, the focus of the 3D forward modeling of CSEMs is reduced to determining the EM fields of ED or MD sources for the background medium. The quasi-minimal residual (QMR) method is used to solve the large sparse complex linear system. Once the FD EM fields have been calculated, the time-domain (TD) response can be obtained using the cosine/sine transformation. The numerical results show that the relative error is less than 5% between the 3D numerical and analytical solutions, which verifies the accuracy of the solver. We further study the difference between the real (bent) and theoretical (straight) wires. We suggest that the shape of the source must be considered for TD and FD CSEMs with a wire source during data processing and inversion. The last example investigated the characteristics of FD EM fields from a finite-length wire and TD EM fields from a rectangular fixed loop on the same conductive tilted disk model buried in resistive sediments. According to the numerical results, we recommend FD CSEMs with a wire source for detecting deep anomalies.

Comparison study of different forward modeling approaches by using the finite-element method for 3D magnetotelluric data in anisotropic media

We have examined three approaches for analyzing 3D magnetotelluric (MT) data with general anisotropy. These included the direct electromagnetic (EM) approach, which solved the electric field Helmholtz equation directly, and the Coulomb-gauged and Ungauged approaches. The three approaches are discretized using finite-element (FE) methods. First, the correctness of the three approaches is validated using a typical isotropic 3D model and a 2D anisotropic model. The FE solutions of each approach obtained by the direct solver MUMPS are compared with the analytical solutions to demonstrate that all approaches generated accurate results. Subsequently, an anisotropic 3D model is designed to investigate the forward characteristics of the different approaches. The FE solutions are obtained using the quasi-minimum residual (QMR) solver with a symmetric successive overrelaxation (SSOR) preconditioner in this model. The forward modeling results of each approach were compared in terms of stiffness matrices, the convergence of the iterative solver, and apparent resistivity and phase. We find that the FE solution of the Coulomb-gauged and Ungauged approaches produces larger stiffness matrices, but the proportion of non-zero elements in the Coulomb-gauged approach is the least. The Ungauged approach can yield faster solutions over a wide range of frequencies. Finally, a complex numerical example is provided to demonstrate the good convergence behavior of the Ungauged approach.

Three-dimensional controlled-source electromagnetic forward modeling by edge-based finite element with a divergence correction

We have developed an edge-based finite-element (FE) modeling algorithm with a divergence correction for calculating the controlled-source electromagnetic (CSEM) responses of a 3D conductivity earth model. We solve a curl-curl equation to directly calculate the secondary electric field to eliminate the source singularity. The choice of the edge-based FE method enables us to properly handle the discontinuity of the normal component of electric fields across conductivity boundaries. Although we can solve the resulting complex-symmetric linear system of equations efficiently by a quasiminimal residual (QMR) method preconditioned with an incomplete Cholesky decomposition for the high-frequency band, the iterative solution process encounters a common problem in the field formulation and does not converge within a practically feasible number of iterations for low frequencies. To overcome this difficulty and to accelerate the iterative solution process in general, we combine a divergence correction technique with the secondary field solution using the QMR solver. We have found that applying the divergence correction intermittently during the iterative solution process ensures the calculation of sufficiently accurate electric and magnetic fields and can significantly speed up the solution process by more than an order of magnitude. We have tested the efficiency and accuracy of our algorithm with 1D and 3D models, and we have found that the divergence correction technique is able to guide the electric field to satisfy the boundary conditions across conductivity interfaces. Although there is a computational overhead required for applying the divergence correction, that cost is significantly offset by the substantial gains in the solution accuracy and speed-up. The work makes the field-based curl-curl formulation using edge elements an efficient and practical method for CSEM simulations.

Research on the forward modeling of controlled-source audio-frequency magnetotelluric in three-dimensional electrical conductivity and magnetic permeability heterogeneity media

We have developed a three-dimensional (3-D) algorithm for computing the controlled-source audio-frequency magnetotelluric response which is capable of handling local inhomogeneities in the electrical conductivity and magnetic permeability distribution. It is a difference equation algorithm that is based on the integral forms of Maxwell’s equations. To avoid a spatial singularity being introduced by the source, we used a secondary electric field approach. We adopt a staggered-grid finite difference method to disperse the secondary electric field equation which is simple to realize and can calculate relatively complex models. The final sparse equation system after dispersion is solved using an ILU preconditioning quasi-minimal residual (QMR) method. The code is validated for nonmagnetic and permeable conductive structures by two synthetic models. One of them is a two-dimensional (2-D) nonmagnetic conductive geoelectric model; we validate our 3-D computation result by comparing it to the result of 2-D CSEM finite element forward program. The other one is a 3-D permeable geoelectric model; we validate our 3-D computation result by using the COMSOL Multiphysics software. We also designed two theoretical 3-D geoelectrical models to analyze the effects of magnetic permeability heterogeneity. Numerical simulation results well reveal the influence of magnetic permeability heterogeneity on high-resistance and low-resistance models.

On the convergence of Q-OR and Q-MR Krylov methods for solving nonsymmetric linear systems

This paper addresses the convergence behavior of Krylov methods for nonsymmetric linear systems which can be classified as quasi-orthogonal (Q-OR) or quasi-minimum residual (Q-MR) methods. It explores, more precisely, whether the influence of eigenvalues is the same when using non-orthonormal bases as it is for the FOM and GMRES methods. It presents parametrizations of the classes of matrices with a given spectrum and right-hand sides generating prescribed Q-OR/Q-MR (quasi) residual norms and discusses non-admissible residual norm sequences. It also gives closed-form expressions of the Q-OR/Q-MR (quasi) residual norms as functions of the eigenvalues and eigenvectors of the matrix of the linear system.

The Non–Symmetric s–Step Lanczos Algorithm: Derivation of Efficient Recurrences and Synchronization–Reducing Variants of BiCG and QMR

Abstract The Lanczos algorithm is among the most frequently used iterative techniques for computing a few dominant eigenvalues of a large sparse non-symmetric matrix. At the same time, it serves as a building block within biconjugate gradient (BiCG) and quasi-minimal residual (QMR) methods for solving large sparse non-symmetric systems of linear equations. It is well known that, when implemented on distributed-memory computers with a huge number of processes, the synchronization time spent on computing dot products increasingly limits the parallel scalability. Therefore, we propose synchronization-reducing variants of the Lanczos, as well as BiCG and QMR methods, in an attempt to mitigate these negative performance effects. These so-called s-step algorithms are based on grouping dot products for joint execution and replacing time-consuming matrix operations by efficient vector recurrences. The purpose of this paper is to provide a rigorous derivation of the recurrences for the s-step Lanczos algorithm, introduce s-step BiCG and QMR variants, and compare the parallel performance of these new s-step versions with previous algorithms.

A new quasi-minimal residual method based on a biconjugate A-orthonormalization procedure and coupled two-term recurrences

Recently, some novel variants of Lanczos-type methods were explored that are based on the biconjugate A-orthonormalization process. Numerical experiments coming from some physical problems indicate that these new methods are competitive with or superior to other popular Krylov subspace methods. Among them, the biconjugate A-orthogonal residual (BiCOR) method is the archetype method. However, like the biconjugate gradient (BiCG) method, the BiCOR method often shows irregular convergence behavior which can lead to numerical instability. To overcome this drawback, motivated by the effectiveness and robustness of the quasi-minimal residual (QMR) method, we derive a new QMR-like approach based on the coupled two-term biconjugate A-orthonormalization process and simple recurrences instead of the QR decomposition. Some convergence properties are given, and the special case for solving complex symmetric linear system is considered. Finally, numerical experiments are reported to illustrate the performances of our methods and their preconditioners on linear systems discretizing the Helmholtz equation or taken from Florida Sparse Matrix Collection.

General constraint preconditioning iteration method for singular saddle-point problems

For the singular saddle-point problems with nonsymmetric positive definite (1,1) block, we present a general constraint preconditioning (GCP) iteration method based on a singular constraint preconditioner. Using the properties of the Moore–Penrose inverse, the convergence properties of the GCP iteration method are studied. In particular, for each of the two different choices of the (1,1) block of the singular constraint preconditioner, a detailed convergence condition is derived by analyzing the spectrum of the iteration matrix. Numerical experiments are used to illustrate the theoretical results and examine the effectiveness of the GCP iteration method. Moreover, the preconditioning effects of the singular constraint preconditioner for restarted generalized minimum residual (GMRES) and quasi-minimal residual (QMR) methods are also tested.

Global Quasi-Minimal Residual Method for Image Restoration

The global quasi-minimal residual (QMR) method is a popular iterative method for the solution of linear systems with multiple right-hand sides. In this paper, we consider the application of the global QMR method to classical ill-posed problems arising from image restoration. Since the scale of the problem is usually very large, the computations with the blurring matrix can be very expensive. In this regard, we use a Kronecker product approximation of the blurring matrix to benefit the computation. In order to reduce the disturbance of noise to the solution, the Tikhonov regularization technique is adopted to produce better approximation of the desired solution. Numerical results show that the global QMR method outperforms the classic CGLS method and the global GMRES method.

Quasi-Minimal Residual Variants of the COCG and COCR Methods for Complex Symmetric Linear Systems in Electromagnetic Simulations

The conjugate orthogonal conjugate gradient (COCG) method has been considered an attractive part of the Lanczos-type Krylov subspace method for solving complex symmetric linear systems. However, it is often faced with apparently irregular convergence behaviors in practical electromagnetic simulations. To avoid such a problem, the symmetric quasi-minimal residual (QMR) method has been developed. On the other hand, the conjugate A-orthogonal conjugate residual (COCR) method, which can be regarded as an extension of the conjugate residual method, also had been established. It shows that the COCR often gives smoother convergence behavior than the COCG method. The purpose of this paper is to apply the QMR approaches to the COCG and COCR to derive two new methods (including their preconditioned versions), and to report the benefits of the modified methods by some practical examples arising in electromagnetic simulations.

A quasi-minimal residual variant of IDRstab using the residual smoothing technique

Quasi-Minimal Residual (QMR)-type methods have been developed for improving the oscillations in the residual norms of the Bi-Conjugate Gradient (Bi-CG) method and the hybrid Bi-CG methods. The IDRstab method, which combines the Induced Dimension Reduction (IDR)(s) method with higher-order stabilizing polynomials, is often more effective than IDR(s) and the hybrid Bi-CG methods, but the residual norms oscillate. In this paper, we design a QMR variant of IDRstab by using a residual smoothing technique in order to obtain smoother convergence behavior. Numerical experiments show the effectiveness of our QMR variant of IDRstab.

Electromagnetic divergence correction for 3D anisotropic EM modeling

Until now, the application of divergence correction has been focusing on the electromagnetic (EM) isotropic modeling. However, in regions where the earth demonstrates strong electrical anisotropy due to lamina or faults saturated with water, the traditional divergence corrections may not be working for EM modeling. In this paper, based on the divergence-free property of the current density, we propose a divergence correction technique to speed up the EM modeling process for a 3D arbitrarily anisotropic earth. The volume current density weighting method is adopted for the discretization of divergence correction, so it can be applied to non-uniform grid model. By enforcing the constraints associated with the divergence-free property, the quasi-minimal residual (QMR) iterations and the total computational time for the EM modeling are largely reduced. Furthermore, we discuss how to divide the QMR process into sections and the influence of the number of iterations and the normalized residual on the solution process. We take a helicopter-borne EM system over a 3D earth with arbitrarily electrical anisotropy as an example to examine the effectiveness of our anisotropic divergence correction technique on the EM modelings.

Project APhiD: A Lorenz-gauged A-[formula omitted] decomposition for parallelized computation of ultra-broadband electromagnetic induction in a fully heterogeneous Earth

An essential element for computational hypothesis testing, data inversion and experiment design for electromagnetic geophysics is a robust forward solver, capable of easily and quickly evaluating the electromagnetic response of arbitrary geologic structure. The usefulness of such a solver hinges on the balance among competing desires like ease of use, speed of forward calculation, scalability to large problems or compute clusters, parsimonious use of memory access, accuracy and by necessity, the ability to faithfully accommodate a broad range of geologic scenarios over extremes in length scale and frequency content. This is indeed a tall order. The present study addresses recent progress toward the development of a forward solver with these properties. Based on the Lorenz-gauged Helmholtz decomposition, a new finite volume solution over Cartesian model domains endowed with complex-valued electrical properties is shown to be stable over the frequency range 10−2–1010Hz and range 10−3–105m in length scale. Benchmark examples are drawn from magnetotellurics, exploration geophysics, geotechnical mapping and laboratory-scale analysis, showing excellent agreement with reference analytic solutions. Computational efficiency is achieved through use of a matrix-free implementation of the quasi-minimum-residual (QMR) iterative solver, which eliminates explicit storage of finite volume matrix elements in favor of “on the fly” computation as needed by the iterative Krylov sequence. Further efficiency is achieved through sparse coupling matrices between the vector and scalar potentials whose non-zero elements arise only in those parts of the model domain where the conductivity gradient is non-zero. Multi-thread parallelization in the QMR solver through OpenMP pragmas is used to reduce the computational cost of its most expensive step: the single matrix–vector product at each iteration. High-level MPI communicators farm independent processes to available compute nodes for simultaneous computation of multi-frequency or multi-transmitter responses.

Efficient pre-conditioned iterative solution strategies for the electromagnetic diffusion in the Earth: finite-element frequency-domain approach

We formulate a 3-D finite-element frequency-domain (FEFD) solution for electromagnetic (EM) diffusion and present efficient solution strategies. A system of FEFD equations is pre-conditioned by incomplete LU (ILU) and subsequently solved by the quasi-minimal residual (QMR) method. A rule of thumb for choosing an effective drop tolerance of ILU is proposed. When multiple sources are simulated in a given model, ILU is computed only once and is reused as a pre-conditioner for multiple QMR computations with different source vectors. Resulting solution vectors are also bootstrapped to reduce the number of QMR iterations required for the convergence. We demonstrate that when conductivity structures of an earth model and source frequencies are updated/perturbed, ILU that is computed from the previous model is still an effective and useful pre-conditioner for new forward modelling problems. Using the reusability of ILU, we also propose a new efficient way to overcome the slow convergence rate of the iterative FEFD solution in the static limit. We show that the reuse of ILU and solution bootstrapping serve as effective strategies for improving the computational efficiency of the iterative FEFD solution. Finally, we apply the proposed efficient solution strategies to marine EM survey scenarios in complex offshore models and further demonstrate their effectiveness.

A new Lanczos-based algorithm for simulating high-frequency two-dimensional electron spin resonance spectra

The Lanczos algorithm (LA) is a useful iterative method for the reduction of a large matrix to tridiagonal form. It is a storage efficient procedure requiring only the preceding two Lanczos vectors to compute the next. The quasi-minimal residual (QMR) method is a powerful method for the solution of linear equation systems, Ax = b. In this report we provide another application of the QMR method: we incorporate QMR into the LA to monitor the convergence of the Lanczos projections in the reduction of large sparse matrices. We demonstrate that the combined approach of the LA and QMR can be utilized efficiently for the orthogonal transformation of large, but sparse, complex, symmetric matrices, such as are encountered in the simulation of slow-motional 1D- and 2D-electron spin resonance (ESR) spectra. Especially in the 2D-ESR simulations, it is essential that we store all of the Lanczos vectors obtained in the course of the LA recursions and maintain their orthogonality. In the LA-QMR application, the QMR weight matrix mitigates the problem that the Lanczos vectors lose orthogonality after many LA projections. This enables substantially more Lanczos projections, as required to achieve convergence for the more challenging ESR simulations. It, therefore, provides better accuracy for the eigenvectors and the eigenvalues of the large sparse matrices originating in 2D-ESR simulations than does the previously employed method, which is a combined approach of the LA and the conjugate-gradient (CG) methods, as evidenced by the quality and convergence of the 2D-ESR simulations. Our results show that very slow-motional 2D-ESR spectra at W-band (95 GHz) can be reliably simulated using the LA-QMR method, whereas the LA-CG consistently fails. The improvements due to the LA-QMR are of critical importance in enabling the simulation of high-frequency 2D-ESR spectra, which are characterized by their very high resolution to molecular orientation.

The Crout version of ILU preconditioning applied to the fast‐multipole method for solving large‐scale electromagnetic open PEC problems

PurposeThe purpose of this paper is to investigate and analyze the efficiency and stability of the implementation of the Crout version of ILU (ILUC) preconditioning on fast‐multipole method (FMM) for solving large‐scale dense complex linear systems arising from electromagnetic open perfect electrical conductor (PEC).Design/methodology/approachThe FMM is employed to reduce the computational complexity of the matrix‐vector product and the memory requirement of the impedance matrix. The numerical examples are initially solved by the quasi‐minimal residual (QMR) method with ILUC preconditioning. In order to fully investigate the performance of ILUC in connection with other iterative solvers, a case is also solved by bi‐conjugate gradient solver and conjugate gradient squared solver with ILUC preconditioning.FindingsThe solutions show that the ILUC preconditioner is stable and significantly improves the performance of the QMR solver on large ill‐conditioned open PEC problems compared to using ILU(0) and threshold‐based ILU (ILUT) preconditioners. It dramatically decreases the number of iterations required for convergence and consequently reduces the total CPU solving time with a reasonable overhead in memory.Practical implicationsThe preconditioning scheme can be applied to large ill‐conditioned open PEC problems to effectively speed up the overall electromagnetic simulation progress while maintaining the computational complexity of FMM. More complex structures including wire‐PEC junctions and microstrip arrays may be addressed in future work.Originality/valueThe performance of ILUC has been previously reported only on preconditioning sparse linear systems, in which the ILU preconditioner is constructed by the ILUC of the coefficient matrix (e.g. matrix arised from two‐dimensional finite element convection‐diffusion problem) and subsequently applied to the same sparse linear systems; so it is important to report its performance on the dense complex linear systems that arised from open PEC electromagnetic problems. In contrast, the preconditioner is constructed upon the near‐field matrix of the FMM and subsequently applied to the whole dense linear system. The comparison of its performance against the diagonal, ILU(0) and ILUT precoditioners is also presented.

Numerical experiments of preconditioned Krylov subspace methods solving the dense non-symmetric systems arising from BEM

Discretization of boundary integral equations leads, in general, to fully populated non-symmetric linear systems of equations. An inherent drawback of boundary element method (BEM) is that, the non-symmetric dense linear systems must be solved. For large-scale problems, the direct methods require expensive computational cost and therefore the iterative methods are perhaps more preferable. This paper studies the comparative performances of preconditioned Krylov subspace solvers as bi-conjugate gradient (Bi-CG), generalized minimal residual (GMRES), conjugate gradient squared (CGS), quasi-minimal residual (QMR) and bi-conjugate gradient stabilized (Bi-CGStab) for the solution of dense non-symmetric systems. Several general preconditioners are also considered and assessed. The results of numerical experiments suggest that the preconditioned Krylov subspace methods are effective approaches solving the large-scale dense non-symmetric linear systems arising from BEM.

2D marine controlled-source electromagnetic modeling: Part 1 — An adaptive finite-element algorithm

In Part 1 of this work, we develop an adaptive finite-element algorithm for forward modeling of the frequency-domain, marine controlled-source electromagnetic (CSEM) response of a 2D conductivity structure that is excited by a horizontal electric dipole source. After transforming the governing equations for the secondary electromagnetic fields into the wavenumber domain, the coupled system of two partial differential equations for the strike-parallel electric and magnetic fields is approximated using the finite-element method. The model domain is discretized using an unstructured triangular element grid that readily accommodates arbitrarily complex structures. A numerical solution of the system of linear equations is obtained using the quasi-minimal residual (QMR) method, which requiresmuch less storagethan full matrix inversion methods. We implement an automatedadaptive grid refinement algorithm in which the finite-element solution is computed iteratively on successively refined grids. Grid refinement is guided by an a posteriori error estimator based on a recently developed gradient recovery operator. The error estimator uses the solution to a dual problem in order to bias refinement toward elements that affect the solution at the electromagnetic (EM) receiver locations and enables the computation of asymptotically exact solutions to the 2.5D partial differential equations. We validate the finite element formulation against the canonical 1D reservoir model and study the performance of the adaptive refinement algorithm. An example model study of a complex offshore structure of interest for petroleum exploration illustrates the utility of the adaptive finite-element method for CSEM modeling.

Global quasi-minimal residual method for the Sylvester equations

In this paper, a global quasi-minimal residual (QMR) method was presented for solving the Sylvester equations. Some properties were investigated with a new matrix product for the global QMR method. Numerical results with the global QMR and GMRES methods compared with the block GMRES method were given. The results show that the global QMR method is less time-consuming than the global GMRES (generalized minimal residual) and block GMRES methods in some cases.

Numerical implementation of the QMR algorithm by using discrete stochastic arithmetic

In each step of the quasi-minimal residual (QMR) method which uses a look-ahead variant of the nonsymmetric Lanczos process to generate basis vectors for the Krylov subspaces induced byA, it is necessary to decide whether to construct the Lanczos vectorsvn+1 andwn+1 as regular or inner vectors. For a regular step it is necessary thatDk=WkTVk is nonsingular. Therefore, in the floating-point arithmetic, the smallest singular value of matrix Dk,σmin(Dk), is computed and an inner step is performed ifσmin(Dk)<∈, where ∈ is a suitably chosen tolerance. In practice it is absolutely impossible to choose correctly the value of the tolerance ∈. The subject of this paper is to show how discrete stochastic arithmetic remedies the problem of this tolerance, as well as the problem of the other tolerances which are needed in the other checks of the QMR method with the estimation of the accuracy of some intermediate results. Numerical examples are used to show the good numerical properties.