Conjugate Gradient Solver Research Articles

Contrast and spatial resolution in MREIT using low amplitude current

Magnetic resonance–electrical impedance tomography employs low amplitude currents injected or induced inside an object. The additional magnetic field due to these currents results in a phase in the MR images. In this study, a modified fast spin–echo sequence was used to measure this magnetic field, which is obtained by scaling the MR phase image. A finite element method with first order triangular elements was used for the solution of the forward problem. An iterated sensitivity matrix-based algorithm was developed for the inverse problem. The resulting ill-conditioned matrix equation was regularized using the Tikhonov method and solved using a conjugate gradient solver. The spatial and contrast resolution of the technique was tested using agarose gel phantoms. A circular phantom with 7 cm diameter and 1 cm thickness is used in the phantom experiments. The amplitude of the injected current was 1 mA. 3, 5 and 8 mm diameter insulators and high conductor objects are used for the spatial resolution study and an average full-width half-maximum value of 4.7 mm is achieved for the 3 mm insulator case. For the contrast analysis, the conductivity of a 15 mm object is varied between 44% and 500% with respect to the background and results are compared to the ideal reconstruction.

Start vector generation for implicit Newmark time integration of the wave equation

Different extrapolation strategies are presented for the construction of starting values in implicit Newmark time-stepping schemes used for the solution of the electromagnetic wave equation. The implicit character of this schemes implies the solution of a system of linear equations in each time step. The choice of a good start value derived from the extrapolation of the equation system is shown to improve the convergence of the iterative conjugate gradient solver which reduces the total simulation time

Solutions, algorithms and inter-relations for local minimization search geophysical inversion

This paper presents in a general form the most popular local minimization search solutions for geophysical inverse problems—the Tikhonov regularization solutions, the smoothest model solutions and the subspace solutions, from which the inter-relationships between these solutions are revealed. For the Tikhonov regularization solution, a variety of forms exist—the general iterative formula, the iterative linearized scheme, the Levenberg–Marquardt version, the conjugate gradient solver (CGS) and local-search quadratic approximation CGS. It is shown here that the first three solutions are equivalent and are just a specified form of the gradient solution. The local-search quadratic approximation CGS is shown to be a more general form from which these three solutions emerge. The smoothest model solution (Occam's inversion) is a subset of the Tikhonov regularization solutions, for which a practical algorithm generally comprises two parts: (1) determine the subset of the Tikhonov regularization solutions which satisfy the desired tolerance of data fit, (2) choose the solution which best fits the data in the subset. The solution procedure is actually an adaptive Tikhonov regularization solution. The subspace solution is mathematically no more than a dimension-decreased transform of model parameterization in the Tikhonov regularization solution or the smoothest model solution. The crucial step is to choose the basis vectors which span the whole model space in the parameterization transform. A simple form of the transform is developed which yields greater flexibility for geophysical inverse applications.

A fast algebraic multigrid preconditioned conjugate gradient solver

This work presents a new approach for selecting the coarse grids allowing a faster algebraic multigrid (AMG) preconditioned conjugate gradient solver. This approach is based on an appropriate choice of the parameter α considering the matrix density during the coarsening process which implies in a significant reduction in the matrix dimension at all AMG levels.

Algorithmic optimizations of a conjugate gradient solver on shared memory architectures

OpenMP is an architecture-independent language for programming in the shared memory model. OpenMP is designed to be simple and powerful in terms of programming abstractions. Unfortunately, the architecture-independent abstractions sometimes come with the price of low parallel performance. This is especially true for applications with an unstructured data access pattern running on distributed shared memory systems (DSM). Here, proper data distribution and algorithmic optimizations play a vital role for performance. In this article, we have investigated ways of improving the performance of an industrial class conjugate gradient (CG) solver, implemented in OpenMP running on two types of shared memory systems. We have evaluated bandwidth minimization, graph partitioning and reformulations of the original algorithm reducing global barriers. By a detailed analysis of barrier time and memory system performance, we found that bandwidth minimization is the most important optimization reducing both L2 misses and remote memory accesses. On a uniform memory system, we get perfect scaling. On a NUMA system, the performance is significantly improved with the algorithmic optimizations leaving the system dependent global reduction operations as a bottleneck.

Optimizing the discrete-dipole approximation for sequences of scatterers with identical shapes but differing sizes or refractive indices

We study scattering of light by small particles with identical shapes but either moderately differing sizes or refractive indices by utilizing the discrete-dipole approximation (DDA). Assuming that accurate DDA solutions are available for either a sequence of sizes or refractive indices, we initialize the iterative conjugate gradient solver for a new size or refractive index by making “educated guesses” of the electric field vectors using classical Lagrange, rational-function, and modified Adams–Bashforth–Moulton extrapolation schemes. In the present pilot study, we assess the initialization schemes for spherical and cubic particles. As compared to the common initialization using the incident electric field, we show that careful extrapolation can significantly reduce the number of iterations. At best, the computing time can decrease by an order of magnitude whereas, typically, the improvement is some tens of percent for sizes comparable to the wavelength. In solving large numbers of single-particle scattering problems, initialization via extrapolation can yield substantial savings in computing time. In particular, the present approach should prove useful when the precise scatterer sizes and refractive indices are unknown, e.g., when interpreting astronomical observations of atmosphereless solar-system objects and experimental measurements.

First Arrival Seismic Tomography (FAST) vs. PStomo_eq applied to crooked line seismic data from the Siljan ring area

When seismic profiles deviate significantly from straight lines, the results from 2D traveltime inversion programs will be in error due to the inherent 3D component present in the data. Thus, it is necessary to use a program that can handle the 3D aspects of the acquisition geometry. This study compares the performance and results from two computer programs for 3D seismic tomography. These algorithms are the package for First Arrival Seismic Tomography (FAST) and a Local Earthquake tomography program, PStomo_eq. Although both codes invert for the velocity field using the conjugate gradient solver LSQR, the common smoothness constraint is handled differently. In addition, the programs do not incorporate the same options for user-specified constraints. These differences in implementation are clearly observed in the inverted velocity fields obtained in this study. Both FAST and PStomo_eq are applied to synthetic and real data sets with crooked line geometry. First arrival traveltimes from seismic data acquired in the Siljan ring impact area are used for the real data set test. The results show that FAST gives smoother models than PStomo_eq. On the real data set PStomo_eq showed a better correlation to the information at hand. Different criteria exist for what is desirable in a model; thus, the choice of which program to use will mostly depend on the particular goals of the study.

비정렬격자계를 사용하는 3차원 유동해석코드 개발 (I) - 수치해석방법 -

A conservative pressure-based finite-volume numerical method has been developed for computing flow and heat transfer by using an unstructured grid system. The method admits arbitrary convex polyhedra. Care is taken in the discretization and solution procedures to avoid formulations that are cell-shape-specific. A collocated variable arrangement formulation is developed, i.e. all dependent variables such as pressure and velocity are stored at cell centers. Gradients required for the evaluation of diffusion fluxes and for second-order-accurate convective operators are found by a novel second-order accurate spatial discretization. Momentum interpolation is used to prevent pressure checkerboarding and the SIMPLE algorithm is used for pressure-velocity coupling. The resulting set of coupled nonlinear algebraic equations is solved by employing a segregated approach, leading to a decoupled set of linear algebraic equations for each dependent variable, with a sparse diagonally dominant coefficient matrix. These equations are solved by an iterative preconditioned conjugate gradient solver which retains the sparsity of the coefficient matrix, thus achieving a very efficient use of computer resources.

A comparison of different radiomagnetotelluric data inversion methods for buried waste sites

Abstract This study deals with two-dimensional (2D) inversions of synthetic and observed radiomagnetotelluric (RMT) data on typical buried conductive waste sites in Europe, and with the practical aspects of different inversion algorithms. In the inversion calculations, we used smoothing and L2-norm stabilizers and compared the results. The resolution of the geometry of the highly conductive waste site, in particular, was investigated. In the inversion with the L2-norm stabilizer, we used the least-squares solution with singular value decomposition (LSSVD) and conjugate gradient (CG), whereas only the conjugate gradient solver was used in the 2D-inversion with the smoothing stabilizer. The inversion results of the synthetic data showed a better resolution of the geometry of the highly conductive waste site when a L2-norm stabilizer was applied in the inversion; in particular, a better detection of the bottom of the waste deposit was achieved. Additional model studies were carried out using synthetic RMT data in order to investigate the 2D inversion of RMT data observed on a 3D structure; these studies showed that the use of TM mode data yields a better resolution of the structure than joint inversion of TE and TM modes. 2D inversions of RMT data on a waste site near Cologne showed that the inversion of the TM mode could resolve well the geometry, especially the bottom of the waste site, if information about the background conductivity structure was available. In this case study, inversion with the L2-norm stabilizer produced a sharper image of the waste site than inversion with the smoothing stabilizer, as indicated also by the inversion study that used synthetic data.

Robust algorithm for random resistor networks using hierarchical domain structure

In this study we discuss methods for solving random resistor networks and similar problems. We discuss the node elimination method and we demonstrate its equivalence to the Gaussian elimination scheme, finding a good elimination order, which makes the method highly efficient. The transfer matrix method is shown to be a special case of the node elimination method with an ordering that is far from optimal. We compare the performance of these exact methods with a state of the art conjugate gradient solver. In general, the node elimination method is the faster method.

Numerical solution of two-phase flow for the advection-dominated and non-linear case

This work presents a highly efficient numerical scheme for solving immiscible, advection-dominated two-phase flow in heterogeneous porous media. The pressure equation is decoupled from the saturation equation using an IMPES approach, while the advective terms are decoupled from the capillary diffusive terms in the saturation equation through sequential operator splitting. The parabolic and hyperbolic equations are approximated in time by implicit and explicit schemes, respectively. Damped Newton linearization is applied to the implicit non-linear diffusive step. Mixed hybrid finite elements are applied to the global pressure equation and to the regularized capillary diffusion term. For both linear systems arising from the approximation procedure, an AMG preconditioned conjugate gradient solver is used. A finite volume scheme with slope limiter is applied to the advective step. Numerical comparison with standard preconditioners demonstrates the reliability of the proposed AMG-preconditioner. Benchmark examples illustrate the robustness of the method.

An efficient algorithm for modelling progressive damage accumulation in disordered materials

AbstractThis paper presents an efficient algorithm for the simulation of progressive fracture in disordered quasi‐brittle materials using discrete lattice networks. The main computational bottleneck involved in modelling the fracture simulations using large discrete lattice networks stems from the fact that a new large set of linear equations needs to be solved every time a lattice bond is broken. Using the present algorithm, the computational complexity of solving the new set of linear equations after breaking a bond reduces to a simple triangular solves (forward elimination and backward substitution) using the already Cholesky factored matrix. This algorithm using the direct sparse solver is faster than the Fourier accelerated iterative solvers such as the preconditioned conjugate gradient (PCG) solvers, and eliminates the critical slowing down associated with the iterative solvers that is especially severe close to the percolation critical points. Numerical results using random resistor networks for modelling the fracture and damage evolution in disordered materials substantiate the efficiency of the present algorithm. In particular, the proposed algorithm is especially advantageous for fracture simulations wherein ensemble averaging of numerical results is necessary to obtain a realistic lattice system response. Copyright © 2005 John Wiley & Sons, Ltd.

Parallel ILU preconditioning and parallel mesh adaptation with load balancing for general domain decompositions for the Navier–Stokes equations

AbstractA complete framework for solving the incompressible Navier–Stokes equations in parallel is presented. An unstructured mesh is decomposed into non‐overlapping subdomains corresponding to the number of processors. Each subdomain is adaptively refined independently based on local Reynolds number estimates. Computational load is balanced by transferring element‐octrees between subdomains.A parallel conjugate gradient solver with ILU preconditioning is achieved by resolving node dependencies based on mesh structure. Each node is sorted by category giving an a priori pivoting suited for parallel solution. The parallel solver has convergence rates comparable to serial solvers with a similar ILU strategy. Copyright © 2005 John Wiley & Sons, Ltd.

A Vectorized Conjugate-Gradient Solver for Sparse Systems of Algebraic Equations

This article describes a preconditioned conjugate-gradient (CG) solver suitable for large sparsely populated systems of algebraic equations. The solver adopts the diagonally scaled scheme that retains the low storage requirement of the basic CG method, while improving its rate of convergence. Compared to other preconditioning schemes, such as the incomplete Cholesky factorization, the diagonally scaled scheme requires minor modifications to the basic CG method and can be implemented easily. Like the basic method, the diagonally scaled scheme is readily vectorizable and can perform better than the incomplete Cholesky factorization on vector processors. The article evaluates the performance of the solver by using the standard test case of a flow in a differentially heated square cavity at different values of the Rayleigh number. It also shows the computer-time requirements of the solver on two vector processors, in both scalar and vectorized modes, which clearly show the benefit of vectorization.

A Vectorized Conjugate-Gradient Solver for Sparse Systems of Algebraic Equations

This article describes a preconditioned conjugate-gradient (CG) solver suitable for large sparsely populated systems of algebraic equations. The solver adopts the diagonally scaled scheme that retains the low storage requirement of the basic CG method, while improving its rate of convergence. Compared to other preconditioning schemes, such as the incomplete Cholesky factorization, the diagonally scaled scheme requires minor modifications to the basic CG method and can be implemented easily. Like the basic method, the diagonally scaled scheme is readily vectorizable and can perform better than the incomplete Cholesky factorization on vector processors. The article evaluates the performance of the solver by using the standard test case of a flow in a differentially heated square cavity at different values of the Rayleigh number. It also shows the computer-time requirements of the solver on two vector processors, in both scalar and vectorized modes, which clearly show the benefit of vectorization.

Optimization of the sources in local hyperthermia using a combined finite element-genetic algorithm method

This article describes an optimization process specially designed for local and regional hyperthermia in order to achieve the desired specific absorption rate in the patient. It is based on a genetic algorithm coupled to a finite element formulation. The optimization method is applied to real human organs meshes assembled from computerized tomography scans. A 3D finite element formulation is used to calculate the electromagnetic field in the patient, achieved by radiofrequency or microwave sources. Space discretization is performed using incomplete first order edge elements. The sparse complex symmetric matrix equation is solved using a conjugate gradient solver with potential projection pre-conditionning. The formulation is validated by comparison of calculated specific absorption rate distributions in a phantom to temperature measurements. A genetic algorithm is used to optimize the specific absorption rate distribution to predict the phases and amplitudes of the sources leading to the best focalization. The objective function is defined as the specific absorption rate ratio in the tumour and healthy tissues. Several constraints, regarding the specific absorption rate in tumour and the total power in the patient, may be prescribed. Results obtained with two types of applicators (waveguides and annular phased array) are presented and show the faculties of the developed optimization process.

Parallel Multigrid Preconditioner for the Cardiac Bidomain Model

The bidomain equations are widely used for the simulation of electrical activity in cardiac tissue but are computationally expensive, limiting the size of the problem which can be modeled. The purpose of this study is to determine more efficient ways to solve the elliptic portion of the bidomain equations, the most computationally expensive part of the computation. Specifically, we assessed the performance of a parallel multigrid (MG) preconditioner for a conjugate gradient solver. We employed an operator splitting technique, dividing the computation in a parabolic equation, an elliptical equation, and a nonlinear system of ordinary differential equations at each time step. The elliptic equation was solved by the preconditioned conjugate gradient method, and the traditional block incomplete LU parallel preconditioner (ILU) was compared to MG. Execution time was minimized for each preconditioner by adjusting the fill-in factor for ILU, and by choosing the optimal number of levels for MG. The parallel implementation was based on the PETSc library and we report results for up to 16 nodes on a distributed cluster, for two and three dimensional simulations. A direct solver was also available to compare results for single processor runs. MG was found to solve the system in one third of the time required by ILU but required about 40% more memory. Thus, MG offered an attractive tradeoff between memory usage and speed, since its performance lay between those of the classic iterative methods (slow and low memory consumption) and direct methods (fast and high memory consumption). Results suggest the MG preconditioner is well suited for quickly and accurately solving the bidomain equations.

A hybrid condensed finite element model with GPU acceleration for interactive 3D soft tissue cutting

AbstractTo meet the requirement of computer‐aided medical operations, apart from the real‐time deformation, it is also necessary in the design to simulate the tissue cutting and suturing in a surgery simulation. In this paper, we present a model on topology change and deformation of soft tissue, referred to as the hybrid condensed finite element model, based on the volumetric finite element method. The most important advantage of our model is its ability to achieve an interactive frame rate for the topology change in surgical simulation on standard PC platform. This is achieved through two innovations. One is to apply the condensation technique, by fully calculating the volumetric deformation in the operation part while only calculating the surface nodes in the non‐operation part. Secondly, the major calculation work in the Conjugate Gradient solver for cutting and deformation is migrated from the CPU to the contemporary GPU to promote the calculation. Test examples have been given to show the feasibility and efficiency of the model. Copyright © 2004 John Wiley & Sons, Ltd.

Implicit purification for temperature-dependent density matrices

An implicit purification scheme is proposed for calculation of the temperature-dependent, grand canonical single-particle density matrix, given as a Fermi-Dirac operator expansion in terms of the Hamiltonian. The computational complexity is shown to scale with the logarithm of the polynomial order of the expansion, or equivalently, with the logarithm of the inverse temperature. The system of linear equations that arise in each implicit purification iteration is solved efficiently by a conjugate gradient solver. The scheme is particularly useful in connection with linear scaling electronic structure theory based on sparse matrix algebra. The efficiency of the implicit temperature expansion technique is analyzed and compared to some explicit purification methods for the zero temperature density matrix.

Anisotropic resistivity inversion

We present an inversion method for 3D electrical imaging in media with aninhomogeneous and anisotropic conductivity distribution. The conductivitydistribution is discretized via finite elements and is described by a second-ordertensor at each finite element node. The inversion method is formulated as afunctional optimization with an error functional containing terms measuringdata misfit and model covariance by means of smoothness, anisotropyand deviation from a starting model. Including the model covarianceinformation overcomes the problem of ill-posedness at the expense oflimiting the allowed models to the class of models which are compatible withthe provided model covariance information. The discretized form of theerror functional is minimized by a Levenberg–Marquardt type methodusing an iterative preconditioned conjugate gradient solver. The use ofan iterative solver allows one to bypass the actual computation of theJacobian or an inverse system matrix. The use of a memory efficient iterativesolver together with the implementation on parallel computers allowslarge-scale inverse problems, comprising several hundred thousand nodes withhundreds of sources and receivers, to be solved. The new method is testedusing computer-generated data from two- and three-dimensional syntheticmodels. For each inversion a choice of penalty parameters, gauging thelevel of model covariance information imposed, has to be made and thelevel of regularization required is hard to estimate. We find that runninga suite of inversions with varying penalty parameters and subsequentexamination of the results (including inspection of residual maps) offers a viablemethod for choosing appropriate numerical values for the penalty levels. Inthe applications we found the inversion process to be highly non-linear.Inversion models from intermediate steps of the iterative inversion showstructure in places that do not exhibit structure in the true model andonly at later iterations do anomalies move to the correct location in themodelling domain. This result indicates that linearized inversions that fail tore-linearize during the inversion process will fail to find meaningful inversionimages. The inversion images achieved using the new method recoverthe important features of the true models, including the approximatemagnitudes of the conductivity anomalies and the magnitudes and directions ofanisotropy anomalies. The inversion images are generally ‘blurred’, that issharp edges are smoothed, and the recovered magnitudes of conductivity,anisotropy and anisotropy direction are generally under-estimated.