A fast elliptic solver for simply connected domains

We consider high order finite difference approximations to the Helmholtz equation in an exterior domain. We include a simplified absorbing boundary condition to approximate the Sommerfeld radiation condition. This yields a large, but sparse, complex system, which is not self-adjoint and not positive definite. We discretize the equation with a compact fourth or sixth order accurate scheme. We solve this large system of linear equations with a Krylov subspace iterative method. Since the method converges slowly, a preconditioner is introduced, which is a Helmholtz equation but with a modified complex wavenumber. This is discretized by a second or fourth order compact scheme. The system is solved by BICGSTAB with multigrid used for the preconditioner. We study, both by Fourier analysis and computations this preconditioned system especially for the effects of high order discretizations.

The focus of this work is the numerical approximation of differential problems involving partial differential equations (PDEs) of elliptic nature, in the context of modelling and simulation in uid dynamics in porous media. The dissertation aims to contribute with an implementation of a multiscale multigrid algorithm, recently introduced in the literature, designed for solving ill-conditioned large linear systems of equations derived from those classes of PDEs, typically associated with Poisson problems of pressure-velocity with boundary conditions typical of ow in porous media. The Poisson problem discussed here is identied from the coupled convectiondiffusion transport system counterpart of PDEs, with dominated convection, and by a linearization by means the use of an operator splitting approach based on the physics of the underlying model. The methodology used for the discretization of the Poisson elliptic problem is mixed hybrid nite elements. The numerical solution of the resulting linear system will be addressed by a multiscale multigrid preconditioned conjugate gradient (PCG) method. We combine both methodologies in order to capture the distinct wavelengths associated with the different wavelengths from the assosiated self-adjoint Poisson operator, strongly inuenced by the heterogeneity of the geological properties of the porous media, in particular to the absolute permeability, which in turn might exhibit very large uctuations of orders of magnitude. Numerical experiments in applications of uid dynamics problems in porous media are presented and discussed for a verication of the results obtained by direct numerical simulations with the multiscale multigrid algorithm under consideration.

When solving large linear systems of equations arising from the discretization of elliptic boundary value problems, a combination of iterative methods and preconditioners based on incomplete LU factorizations is frequently used. Given a model problem with variable coefficients, we investigate a class of incomplete LU factorizations depending on a relaxation parameter. We show that the associated preconditioner and the factorization itself both are numerically stable. The theoretical results are complemented by numerical experiments.

Finite-Difference Frequency-Domain methods (FDFD) require solution of large linear systems of equations. These large systems are represented by matrix equations including highly sparse coefficient matrices, and they can often only be solved by using iterative methods. This paper presents an algorithm in which the matrix-equation solution approach in an iterative method is replaced by a multi-step solution process. Instead of using a coefficient matrix, the coefficients in the FDFD formulations are kept as three-dimensional arrays, and they are treated as operators. The algorithm is used together with the Bi-Conjugate Gradients Stabilized (BICGSTAB) method. This is applied to a three-dimensional FDFD method to solve for scattering from dielectric objects. It is also applied to two other FDFD methods (a single-grid and a double-grid FDFD) to solve for scattering from chiral objects. It has been shown that the presented algorithm effectively reduces the solution time and memory requirements.

In this paper finite element method for 3D DC resistivity modeling accelerated using multi-GPU (Graphics Processing Unit). Solution of the large system of linear equations is the most expensive computation in finite element method performed in GPUs to reduce the computational time. Conjugate gradient solver used to solve large system of linear equations. We developed kernel for conjugate gradient solver that exploit data vectorization and written in PTX assembly form. We perform test on GTX 750Ti GPU and Tesla C2050 GPU. Our kernel have better sparse matrix-vector performance than CUSPARSE library in the first GPU, but lower performance in the latter GPU. The performance comparisons to the library are about 1.4 times and 0.7 times respectively. Our multi-GPU implementation achieved about 1.9 times performance of single GPU by using 2 identical GPUs. In comparison to the serial CPU implementation, about 10 times speedup could be achieved by using 2 GPUs.

Randomized scalable vector algorithms for calculation of matrix iterations and solving extremely large linear algebraic equations are developed. Among applications presented in this paper are randomized iterative methods for large linear systems of algebraic equations governed by M-matrices. The crucial idea of the randomized method is that the iterations are performed by sampling random columns only, thus avoiding not only matrix-matrix but also matrix-vector multiplications. The suggested vector randomized methods are highly efficient for solving linear equations of high dimension, the computational cost depends only linearly on the dimension.

The numerical simulation of shale gas reservoirs requires transforming mathematical models into large linear equation systems, and solving large linear equation systems requires a lot of time. A novel ILU-GMRES method for solving large linear equations is presented, which treats the sub matrix discretized at nodes as the basic elements of the coefficient matrix. By combining the incomplete LU decomposition method (ILU) with the generalized minimum residue method (GMRES), the discretized large linear equations are iteratively solved. A comparison was made between the new method and other methods for calculating the actual shale gas reservoir model. The results reveal that the new method has a faster calculation speed in solving linear equation systems. The calculation time of the novel ILU-GMRES method is 1029s, the calculation time of the generalized minimum residual method (GMRES) method is 2709s, the calculation time of the Gauss–Seidel method (GS) is 3786s, and the calculation time of the successive over relaxation method (SOR) method is 3386s, the calculation time of the preprocessing conjugate gradient method (PCG) is 3115s. The calculation speed of the novel method is almost 2.6 times faster than the GMRES method, the novel method is more effective in the numerical simulation of shale gas.

Use of the stochastic Galerkin finite element methods leads to large systems\nof linear equations obtained by the discretization of tensor product solution\nspaces along their spatial and stochastic dimensions. These systems are\ntypically solved iteratively by a Krylov subspace method. We propose a\npreconditioner which takes an advantage of the recursive hierarchy in the\nstructure of the global matrices. In particular, the matrices posses a\nrecursive hierarchical two-by-two structure, with one of the submatrices block\ndiagonal. Each one of the diagonal blocks in this submatrix is closely related\nto the deterministic mean-value problem, and the action of its inverse is in\nthe implementation approximated by inner loops of Krylov iterations. Thus our\nhierarchical Schur complement preconditioner combines, on each level in the\napproximation of the hierarchical structure of the global matrix, the idea of\nSchur complement with loops for a number of mutually independent inner Krylov\niterations, and several matrix-vector multiplications for the off-diagonal\nblocks. Neither the global matrix, nor the matrix of the preconditioner need to\nbe formed explicitly. The ingredients include only the number of stiffness\nmatrices from the truncated Karhunen-Lo\\`{e}ve expansion and a good\npreconditioned for the mean-value deterministic problem. We provide a condition\nnumber bound for a model elliptic problem and the performance of the method is\nillustrated by numerical experiments.\n

We develop an analytical and numerical methodology for the analysis of large bottom-hole temperature (BHT) data sets from sedimentary basins, and test the methodology using temperature, stratigraphic, and lithologic data from 411 boreholes in the Michigan Basin. Least-squares estimates of temperature gradients in the formations and lithologies present are calculated as solutions to a large system of linear equations. At each borehole the temperature difference between the bottom and top of the hole is represented as a sum of temperature increments through the various formations or lithologies penetrated by the borehole. Quadratic programming techniques enable bounds to be placed on the gradient solutions in order to suppress or exclude improbable gradient estimates. Numerical experiments with synthetic data reveal that the estimates of temperature gradients for a given formation or lithology are sensitive to the degree of representation of that unit; well represented units have more stable gradient estimates in the presence of noise than do poorly represented units. The estimates of temperature gradients obtained for lithologies are more stable than those for formations and are believed to be good estimates of actual lithologic temperature gradients in the Michigan Basin. Formation temperature gradients obtained as a weighted sum of the estimates of the component lithologic temperature gradients appear to be good estimates of the average temperature gradients for the formations of the basin.At each borehole a temperature residual exists corresponding to the difference between the observed BHT and the BHT predicted by the estimated interval temperature gradients. Residuals are far more stable than estimated temperature gradients. The values of residuals change little regardless of whether lithology, formation, bounded, or unbounded gradient estimates are used to calculate them. Maps of residuals indicate well-defined and spatially coherent patterns of positive and negative temperature residuals. Filtered subsets of large-magnitude residuals alone show a pattern of negative residuals coinciding with the mid-Michigan gravity high, a geophysical feature thought to delineate a Precambrian (Keweenawan) rift zone in the crust beneath the basin. Thermal models of the Michigan Basin and the crust and upper mantle beneath the basin indicate that the suspected rift beneath the basin can cause a variation in basement heat flow sufficient to produce temperature residuals of the magnitude observed in the sediments, with negative temperature residuals over the area of the rift.

Many problems in practice require the solution of very large systems of linear equations Ax = b in which the matrix A, fortunately, is sparse, i.e., has relatively few nonvanishing elements. Systems of this type arise, e.g., in the application of difference methods or finite-element methods to the approximate solution of boundary-value problems in partial differential equations. The usual elimination methods (see Chapter 4) cannot normally be applied here, since without special precautions they tend to lead to the formation of more or less dense intermediate matrices, making the number of arithmetic operations necessary for the solution much too large, even for present-day computers, not to speak of the fact that the intermediate matrices no longer fit into the usually available computer memory.

Iterative Tikhonov regularization of tensor equations based on the Arnoldi process and some of its generalizations

On the application of orthogonal polynomials to the iterative solution of linear systems of equations with indefinite or non-Hermitian matrices

In most photoacoustic tomography (PAT) reconstruction approaches, it is assumed that the receiving transducers have omnidirectional response and can fully surround the region of interest. These assumptions are not satisfied in practice. To deal with these limitations, we present a novel deconvolution based photoacoustic reconstruction with sparsity regularization (DPARS) technique. The DPARS algorithm is a semi-analytical reconstruction approach in which the projections of the absorber distribution derived from a deconvolution-based method are computed and used to generate a large linear system of equations. In these projections, computed over limited viewing angles, the directivity effect of the transducer is taken into account. The distribution of absorbers is computed using a sparse representation of absorber coefficients obtained from the discrete cosine transform. This sparse representation helps improve the numerical conditioning of the system of equations and reduces the computation time of the deconvolution-based approach by one order of magnitude relative to Tikhonov regularization. The algorithm has been tested in simulations, and using two-dimensional and three-dimensional experimental data obtained with a conventional ultrasound transducer. The results show that DPARS, when evaluated using contrast-to-noise ratio and root-mean-square errors, outperforms the conventional delay-and-sum (DAS) reconstruction method.

Sparsified Randomization Monte Carlo (SRMC) algorithms for solving systems of linear algebraic equations introduced in our previous paper [34] are discussed here in a broader context. In particular, I present new randomized solvers for large systems of linear equations, randomized singular value (SVD) decomposition for large matrices and their use for solving inverse problems, and stochastic simulation of fields. Stochastic projection methods, which I call here random row action algorithms, are extended to problems which involve systems of equations and constrains in the form of systems of linear inequalities.

The problem of finding a sparse solution to large systems of linear equations arises in many applications related to signal processing. Sometimes, the size of these systems is so large that the known methods are inefficient. Such systems can be solved only if there is additional structure inherent in them. In this paper, an efficient approach for finding sparse solutions to super-large systems of linear equations with a tensor structure of a certain type is proposed. The theoretical analysis and experimental results make it possible to judge the efficiency of the proposed method.