Quasi-minimal Residual Method Research Articles

Transpose-free quasi-minimal residual method based on tensor format for generalized coupled Sylvester tensor equations

Transpose-free quasi-minimal residual method based on tensor format for generalized coupled Sylvester tensor equations

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.

Fast iterative solver for 3D finite-element modeling of controlled-source electromagnetic responses with vector and scalar potentials

The controlled-source electromagnetic (CSEM) forward modeling with vector and scalar potentials could provide insights into the inductive and galvanic effects. We have developed a finite-element (FE) forward modeling algorithm for the CSEM with vector and scalar potentials. The first-order vector and nodal shape functions are used for the vector and scalar potentials, respectively. To mitigate the null space and thus the nonuniqueness of the potentials caused by the curl operator, we introduce a new preconditioner constructed from an incomplete Cholesky decomposition with zero fill-ins of the Laplacian approximation to the resulted matrix after the FE discretization. We implemented a preconditioned quasiminimal residual method to iteratively solve the resulting linear system with the new preconditioner. We first verified the accuracy of the algorithm through comparison against the analytic solutions obtained for a three-layered earth model. The new iterative modeling algorithm converges much faster compared with the modeling algorithm with the preconditioner constructed from the original stiffness matrix. The efficiency and accuracy are further illustrated by comparison with the modeling scheme based on coupled potentials with explicit enforcement of the Coulomb gauge condition for the vector potential by modeling three 3D models. In addition, the accuracy of our algorithm is demonstrated via a comparison of our numerical solution for a 3D model with the numerical solution obtained from the integral equation method. We also present the inductive and galvanic components of the horizontal secondary electric field for numerical solutions obtained with direct and iterative solvers. The numerical test demonstrated that, aside from the efficiency improved with the new iterative solver, the unstable and nonuniqueness problem for the potentials is eliminated by the new preconditioner with the Laplacian approximation involved.

Recycling Newton–Krylov algorithm for efficient solution of large scale power systems

Power flow calculations are crucial for the study of power systems, as they can be used to calculate bus voltage magnitudes and phase angles, as well as active and reactive power flows on lines. In this paper, a new approach, the Recycling Newton–Krylov (ReNK) algorithm, is proposed to solve the linear systems of equations in Newton–Raphson iterations. The proposed method uses the Generalized Conjugate Residuals with inner orthogonalization and deflated restarting (GCRO-DR) method within the Newton–Raphson algorithm and reuses the Krylov subspace information generated in previous Newton runs. We evaluate the performance of the proposed method over the traditional direct solver (LU) and iterative solvers (Generalized Minimal Residual Method (GMRES), the Biconjugate Gradient Stabilized Method (Bi-CGSTAB) and Quasi-Minimal Residual Method (QMR)) as the inner linear solver of the Newton–Raphson method. We use different test systems with a number of busses ranging from 300 to 70000 and compare the number of iterations of the inner linear solver (for iterative solvers) and the CPU times (for both direct and iterative solvers). We also test the performance of the ReNK algorithm for contingency analysis and for different load conditions to simulate optimization problems and observe possible performance gains.

An Effective Algorithm for 2D Marine CSEM Modeling in Anisotropic Media Using a Wavelet Galerkin Method

The marine controlled-source electromagnetic method (MCSEM) has attracted considerable attention as an approach to explore marine oil and gas resources and geological structures. This study presents a new wavelet Galerkin method (WGM) to solve the forward modeling problem of 2D MCSEM data incorporating conductivity anisotropy. The method uses Daubechies wavelets that may be differentiated based on the need to solve the governing field equations of MCSEM. A quasi-minimal residual method was adopted by combining an incomplete LU preconditioner to solve the WGM equations. The numerical results were compared with the analytical solution and those obtained by the finite difference and element methods. The results show that the proposed WGM is superior to the finite element and difference methods in terms of computing time and memory requirements. This algorithm can be applied to solve the forward modeling problem of MCSEM. The conductivity anisotropy of the background medium affects the MCSEM response more than the reservoir anisotropy. The match between the modeled results and measured data for the simplified real model demonstrates the necessity for using the anisotropic model to interpret data. Although this study used the proposed algorithm for 2D models, it may also be used for 3D models.

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.

A 3-D Hybrid Cell Boundary Element Method for Time-Harmonic Eddy Current Problems on Multiply Connected Domains

A novel 3-D hybrid formulation for time-harmonic eddy current problems in multiply connected domains is presented. The interior problem (in conductive regions) is discretized by the cell method (CM) in terms of magnetic vector potentials a, whereas the exterior problem (in the unbounded air domain) is discretized by the boundary element method (BEM) in terms of reduced scalar potentials ${\varphi }_{r}$ . Novel topological constraints are derived from magnetostatic energy conservation by using a decomposition of the magnetic field which minimizes the support and the number of cohomology generators. A fast algebraic procedure to pre-compute source fields to handle efficiently current-driven coils is also proposed. The final matrix system can be solved in a limited number of iterations by transpose-free quasi-minimal residual method with symmetric successive over-relaxation preconditioning. Convergence tests show that numerical results are in a very good agreement with third-order finite element method (FEM) on a 2-D axisymmetric model. The CM-BEM with piecewise uniform approximation shows also to be very effective when analyzing fully 3-D test cases, since second-order FEM accuracy is attained even with coarse mesh refinements, by using, however, a much lower number of degrees of freedom compared to FEM.

A new method for computation of eigenvector derivatives with distinct and repeated eigenvalues in structural dynamic analysis

Abstract Determining eigenvector derivatives is a challenging task due to the singularity of the coefficient matrices of the governing equations, especially for those structural dynamic systems with repeated eigenvalues. An effective strategy is proposed to construct a non-singular coefficient matrix, which can be directly used to obtain the eigenvector derivatives with distinct and repeated eigenvalues. This approach also has an advantage that only requires eigenvalues and eigenvectors of interest, without solving the particular solutions of eigenvector derivatives. The Symmetric Quasi-Minimal Residual (SQMR) method is then adopted to solve the governing equations, only the existing factored (shifted) stiffness matrix from an iterative eigensolution such as the subspace iteration method or the Lanczos algorithm is utilized. The present method can deal with both cases of simple and repeated eigenvalues in a unified manner. Three numerical examples are given to illustrate the accuracy and validity of the proposed algorithm. Highly accurate approximations to the eigenvector derivatives are obtained within a few iteration steps, making a significant reduction of the computational effort. This method can be incorporated into a coupled eigensolver/derivative software module. In particular, it is applicable for finite element models with large sparse matrices.

Meshfree Method in Geophysical Electromagnetic Prospecting: The 2D Magnetotelluric Example

As an important supplement and development of traditional methods, the meshfree method has received a great deal of attention in the field of engineering calculation, and has been successfully used to solve many problems which traditional methods have difficulty in solving. However, the application of meshfree method is relatively less in the area of geophysics. In this paper, we apply the meshfree method to the numerical simulation of geophysical electromagnetic prospecting, taking the 2D magnetotelluric as an example and deduce the corresponding meshfree radial point interpolation method (RPIM) equivalent linear equations in detail. The high-efficiency and accurate solutions of large-scale sparse linear equations are solved by the quasi-minimal residual method based on Krylov subspace. The optimal values of the shape parameters are given by numerical experiments. The correctness of the meshfree method is verified by a layered model. The root mean square error of the calculation results is no more than 0.35%, its accuracy is superior to the finite element method. We also compare the meshfree solution with FEM solution by calculating an inclined vein body model, and the calculation results are in good agreement. A continuously changing fault model and undulating terrain model which traditional methods have difficulty in simulating are respectively calculated, the sectional profiles of the apparent resistivity accurately reflect the trend of the anomalies. The meshfree method does not require the complicated mesh generation, and the physical parameters are loaded at a series of points, thus it is especially suitable for the calculation of the complex geological models. With the rapid development of computational science, the meshfree techniques will certainly become a new robust numerical simulation method in geophysical electromagnetic prospecting.

A time-integration method for stable simulation of extremely deformable hyperelastic objects

This paper presents a time integration method for realtime simulation of extremely deformable objects subject to geometrically nonlinear hyperelasticity. In the presented method, the equation of motion of the system is discretized by the backward Euler method, and linearly approximated through the first-order Taylor expansion. The approximate linear equation is solved with the quasi-minimal residual method (QMR), which is an iterative linear equation solver for non-symmetric or indefinite matrices. The solution is then corrected considering the nonlinear term that is omitted at the Taylor expansion. The method does not demand the constitutive law to guarantee the positive definiteness of the stiffness matrix. Experimental results show that the presented method realizes stable behavior of the simulated model under such deformation that the tetrahedral elements are almost flattened. It is also shown that QMR outperforms the biconjugate gradient stabilized method (BiCGStab) in this application.

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.

Study on performance of Krylov subspace methods at solving large-scale contact problems

This paper discusses the performance of preconditioned Krylov subspace methods at solving large-scale finite-element-formulated contact problems based on a finite element analysis program (FEAP) and portable, extensible toolkit for scientific computation (PETSc). Different combinations of the common preconditioners and Krylov subspace methods, namely conjugate gradient, generalized minimal residual, conjugate residual, biconjugate gradient, conjugate gradient square, biconjugate gradient stabilized, and transpose-free quasi-minimal residual method, are tested on the problems with rising size. Numerical experiments are carried out in respect of time consumption, parallelism, memory usage, and stability. The results show time consumption for most of combinations increases linearly with the growth of problem size, and block Jacobi and Jacobi-preconditioned conjugate gradient and conjugate residual methods have better potential for large-scale problems in the given and similar case.

A parallel VOF IB pressure-correction method for simulation of multiphase flows

This paper develops an improved pressure correction method, coupled with volume-of-fluid (VOF) and immersed boundary (IB) methods, to accelerate multiphase flow computations. A momentum interpolation method (MIM) is used to avoid the issue of a checkerboard pressure field on a non-staggered grid system. A transpose-free quasi-minimal residual (TFQMR) method is applied to reduce the CPU time when solving the Poisson equation for the pressure correction procedure. A suitable initial guess for the TFQMR method is investigated to enhance the stability of the matrix solver. The IB method is applied to treat the solid-wall boundary conditions while OpenMP is used to reduce CPU time. Two physical problems – dam breakage and a water wave impacting on a tall structure – are simulated and compared with known data. The improved numerical scheme has CPU speeds at least 8 times faster than the initial numerical scheme for simulating two-phase flows. With regard to the parallel computational efficiency, using six cores with the OpenMP allows the improved scheme to have 77% parallel efficiency. These modifications greatly enhance overall computational accuracy and efficiency.

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.

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 global transpose-free method with quasi-minimal residual strategy for the Sylvester equations

The Sylvester equation has numerous applications in control and system theory, power system, linear algebra, model reduction and signal processing. In the present paper, a global transpose-free quasi-minimal residual method for the Sylvester equations is proposed. The resulting method, with short-term recurrences, does not involve the multiplication with the transpose of the coefficient matrices, and often exhibits smoother convergence behavior than some other existing global methods. Finally, several numerical examples are reported by comparing the proposed method with some global methods.