Dennis C. Smolarski

Implementation of an adaptive algorithm for Richardson's method

Abstract We discuss the implementation of an adaptive algorithm proposed by one of us. The algorithm is a hybrid of the gmres method and Richardsons method. Richardsons method ( RM ) depends on a set of parameters that are computed by minimizing the L 2 norm of a polynomial over the convex hull of eigenvalues. Execution of GMRES yields not only an approximate solution but also the approximate convex hull. RM is used to avoid storing and working with a large number of vectors as GMRES often requires. This method is also advantageous for the solution of large problems. We consider several test problems and compare our algorithm primarily with the conjugate-gradient-squared algorithm, but also with GMRES and to CG (applied to the normal equations). For many (test) problems our algorithm takes roughly 50 percent more work than the conjugate-gradient-squared algorithm, although if the matrix is either preconditioned or indefinite, our algorithm is more efficient. However, our algorithm currently imposes an undesirable burden on the user, who is invited to consider a variety of numerical parameters to manipulate, such as the number of steps of rm , in order to enhance performance: the values we suggest are only empirical.

Why Gaussian quadrature in the complex plane

This paper synthesizes formally orthogonal polynomials, Gaussian quadrature in the complex plane and the bi-conjugate gradient method together with an application. Classical Gaussian quadrature approximates an integral over (a region of) the real line. We present an extension of Gaussian quadrature over an arc in the complex plane, which we call complex Gaussian quadrature. Since there has not been any particular interest in the numerical evaluation of integrals over the long history of complex function theory, complex Gaussian quadrature is in need of motivation. Gaussian quadrature in the complex plane yields approximations of certain sums connected with the bi-conjugate gradient method. The scattering amplitude cTA−1b is an example where A is a discretization of a differential–integral operator corresponding to the scattering problem and b and c are given vectors. The usual method to estimate this is to use cTx(k). A result of Warnick is that this is identically equal to the complex Gaussian quadrature estimate of 1/λ. Complex Gaussian quadrature thereby replaces this particular inner product in the estimate of the scattering amplitude.

On the performance of SPAI and ADI-like preconditioners for core collapse supernova simulations in one spatial dimension

The simulation of core collapse supernovae calls for the time accurate solution of the (Euler) equations for inviscid hydrodynamics coupled with the equations for neutrino transport. The time evolution is carried out by evolving the Euler equations explicitly and the neutrino transport equations implicitly. Neutrino transport is modeled by the multi-group Boltzmann transport (MGBT) and the multi-group flux limited diffusion (MGFLD) equations. An implicit time stepping scheme for the MGBT and MGFLD equations yields Jacobian systems that necessitate scaling and preconditioning. Two types of preconditioners, namely, a sparse approximate inverse (SPAI) preconditioner and a preconditioner based on the alternating direction implicit iteration (ADI-like) have been found to be effective for the MGFLD and MGBT formulations. This paper compares these two preconditioners. The ADI-like preconditioner performs well with both MGBT and MGFLD systems. For the MGBT system tested, the SPAI preconditioner did not give competitive results. However, since the MGBT system in our experiments had a high condition number before scaling and since we used a sequential platform, care must be taken in evaluating these results.

Scalable, Hydrodynamic and Radiation-Hydrodynamic Studies of Neutron Stars Mergers

We discuss the high performance computing issues involved in the numerical simulation of binary neutron star mergers and supernovae. These phenomena, which are of great interest to astronomers and physicists, can only be described by modeling the gravitational field of the objects along with the flow of matter and radiation in a self consistent manner. In turn, such models require the solution of the gravitational field equations, Eulerian hydrodynamic equations, and radiation transport equations. This necessitates the use of scalable, high performance computing assets to conduct the simulations. We discuss some of the parallel computing aspects of this challenging task in this paper.