An iterative semi-implicit scheme with robust damping

An efficient, iterative semi-implicit (SI) numerical method for the time integration of stiff wave systems is presented. Physics-based assumptions are used to derive a convergent iterative formulation of the SI scheme which enables the monitoring and control of the error introduced by the SI operator. This iteration essentially turns a semi-implicit method into a fully implicit method. Accuracy, rather than stability, determines the timestep. The scheme is second-order accurate and shown to be equivalent to a simple preconditioning method. We show how the diffusion operators can be handled so as to yield the property of robust damping, i.e., dissipating the solution at all values of the parameter \mathcal D\dt, where $\mathcal D$ is a diffusion operator and \dt the timestep. The overall scheme remains second-order accurate even if the advection and diffusion operators do not commute. In the limit of no physical dissipation, and for a linear test wave problem, the method is shown to be symplectic. The method is tested on the problem of Kinetic Alfv\'en wave mediated magnetic reconnection. A Fourier (pseudo-spectral) representation is used. A 2-field gyrofluid model is used and an efficacious k-space SI operator for this problem is demonstrated. CPU speed-up factors over a CFL-limited explicit algorithm ranging from $\sim20$ to several hundreds are obtained, while accurately capturing the results of an explicit integration. Possible extension of these results to a real-space (grid) discretization is discussed.Comment: Submitted to the Journal of Computational Physics. Clarifications and caveats in response to referees, numerical demonstration of convergence rate, generalized symplectic proo

Spectral Methods for Numerical Relativity. The Initial Data Problem

Numerical relativity has traditionally been pursued via finite differencing. Here we explore pseudospectral collocation (PSC) as an alternative to finite differencing, focusing particularly on the solution of the Hamiltonian constraint (an elliptic partial differential equation) for a black hole spacetime with angular momentum and for a black hole spacetime superposed with gravitational radiation. In PSC, an approximate solution, generally expressed as a sum over a set of orthogonal basis functions (e.g., Chebyshev polynomials), is substituted into the exact system of equations and the residual minimized. For systems with analytic solutions the approximate solutions converge upon the exact solution exponentially as the number of basis functions is increased. Consequently, PSC has a high computational efficiency: for solutions of even modest accuracy we find that PSC is substantially more efficient, as measured by either execution time or memory required, than finite differencing; furthermore, these savings increase rapidly with increasing accuracy. The solution provided by PSC is an analytic function given everywhere; consequently, no interpolation operators need to be defined to determine the function values at intermediate points and no special arrangements need to be made to evaluate the solution or its derivatives on the boundaries. Since the practice of numerical relativity by finite differencing has been, and continues to be, hampered by both high computational resource demands and the difficulty of formulating acceptable finite difference alternatives to the analytic boundary conditions, PSC should be further pursued as an alternative way of formulating the computational problem of finding numerical solutions to the field equations of general relativity.Comment: 15 pages, 5 figures, revtex, submitted to PR

Composing Scalable Nonlinear Algebraic Solvers

Most efficient linear solvers use composable algorithmic components, with the most common model being the combination of a Krylov accelerator and one or more preconditioners. A similar set of concepts may be used for nonlinear algebraic systems, where nonlinear composition of different nonlinear solvers may significantly improve the time to solution. We describe the basic concepts of nonlinear composition and preconditioning and present a number of solvers applicable to nonlinear partial differential equations. We have developed a software framework in order to easily explore the possible combinations of solvers. We show that the performance gains from using composed solvers can be substantial compared with gains from standard Newton-Krylov methods.Comment: 29 pages, 14 figures, 13 table

A Robust Solution Procedure for Hyperelastic Solids with Large Boundary Deformation

Compressible Mooney-Rivlin theory has been used to model hyperelastic solids, such as rubber and porous polymers, and more recently for the modeling of soft tissues for biomedical tissues, undergoing large elastic deformations. We propose a solution procedure for Lagrangian finite element discretization of a static nonlinear compressible Mooney-Rivlin hyperelastic solid. We consider the case in which the boundary condition is a large prescribed deformation, so that mesh tangling becomes an obstacle for straightforward algorithms. Our solution procedure involves a largely geometric procedure to untangle the mesh: solution of a sequence of linear systems to obtain initial guesses for interior nodal positions for which no element is inverted. After the mesh is untangled, we take Newton iterations to converge to a mechanical equilibrium. The Newton iterations are safeguarded by a line search similar to one used in optimization. Our computational results indicate that the algorithm is up to 70 times faster than a straightforward Newton continuation procedure and is also more robust (i.e., able to tolerate much larger deformations). For a few extremely large deformations, the deformed mesh could only be computed through the use of an expensive Newton continuation method while using a tight convergence tolerance and taking very small steps.Comment: Revision of earlier version of paper. Submitted for publication in Engineering with Computers on 9 September 2010. Accepted for publication on 20 May 2011. Published online 11 June 2011. The final publication is available at http://www.springerlink.co

Eigenvector Model Descriptors for Solving an Inverse Problem of Helmholtz Equation: Extended Materials

We study the seismic inverse problem for the recovery of subsurface properties in acoustic media. In order to reduce the ill-posedness of the problem, the heterogeneous wave speed parameter to be recovered is represented using a limited number of coefficients associated with a basis of eigenvectors of a diffusion equation, following the regularization by discretization approach. We compare several choices for the diffusion coefficient in the partial differential equations, which are extracted from the field of image processing. We first investigate their efficiency for image decomposition (accuracy of the representation with respect to the number of variables and denoising). Next, we implement the method in the quantitative reconstruction procedure for seismic imaging, following the Full Waveform Inversion method, where the difficulty resides in that the basis is defined from an initial model where none of the actual structures is known. In particular, we demonstrate that the method is efficient for the challenging reconstruction of media with salt-domes. We employ the method in two and three-dimensional experiments and show that the eigenvector representation compensates for the lack of low frequency information, it eventually serves us to extract guidelines for the implementation of the method.Comment: 45 pages, 37 figure

Seismic Ray Impedance Inversion

This thesis investigates a prestack seismic inversion scheme implemented in the ray parameter domain. Conventionally, most prestack seismic inversion methods are performed in the incidence angle domain. However, inversion using the concept of ray impedance, as it honours ray path variation following the elastic parameter variation according to Snell’s law, shows the capacity to discriminate different lithologies if compared to conventional elastic impedance inversion. The procedure starts with data transformation into the ray-parameter domain and then implements the ray impedance inversion along constant ray-parameter profiles. With different constant-ray-parameter profiles, mixed-phase wavelets are initially estimated based on the high-order statistics of the data and further refined after a proper well-to-seismic tie. With the estimated wavelets ready, a Cauchy inversion method is used to invert for seismic reflectivity sequences, aiming at recovering seismic reflectivity sequences for blocky impedance inversion. The impedance inversion from reflectivity sequences adopts a standard generalised linear inversion scheme, whose results are utilised to identify rock properties and facilitate quantitative interpretation. It has also been demonstrated that we can further invert elastic parameters from ray impedance values, without eliminating an extra density term or introducing a Gardner’s relation to absorb this term. Ray impedance inversion is extended to P-S converted waves by introducing the definition of converted-wave ray impedance. This quantity shows some advantages in connecting prestack converted wave data with well logs, if compared with the shearwave elastic impedance derived from the Aki and Richards approximation to the Zoeppritz equations. An analysis of P-P and P-S wave data under the framework of ray impedance is conducted through a real multicomponent dataset, which can reduce the uncertainty in lithology identification.Inversion is the key method in generating those examples throughout the entire thesis as we believe it can render robust solutions to geophysical problems. Apart from the reflectivity sequence, ray impedance and elastic parameter inversion mentioned above, inversion methods are also adopted in transforming the prestack data from the offset domain to the ray-parameter domain, mixed-phase wavelet estimation, as well as the registration of P-P and P-S waves for the joint analysis. The ray impedance inversion methods are successfully applied to different types of datasets. In each individual step to achieving the ray impedance inversion, advantages, disadvantages as well as limitations of the algorithms adopted are detailed. As a conclusion, the ray impedance related analyses demonstrated in this thesis are highly competent compared with the classical elastic impedance methods and the author would like to recommend it for a wider application

Solution of the inverse scattering problem by T-matrix completion. II. Simulations

This is Part II of the paper series on data-compatible T-matrix completion (DCTMC), which is a method for solving nonlinear inverse problems. Part I of the series contains theory and here we present simulations for inverse scattering of scalar waves. The underlying mathematical model is the scalar wave equation and the object function that is reconstructed is the medium susceptibility. The simulations are relevant to ultrasound tomographic imaging and seismic tomography. It is shown that DCTMC is a viable method for solving strongly nonlinear inverse problems with large data sets. It provides not only the overall shape of the object but the quantitative contrast, which can correspond, for instance, to the variable speed of sound in the imaged medium.Comment: This is Part II of a paper series. Part I contains theory and is available at arXiv:1401.3319 [math-ph]. Accepted in this form to Phys. Rev.