SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

All Journal content published prior to 1997 is part of LOCUS.

Search Issue | RSS Feeds RSS
Next Issue

2002

Volume 24, Issue 1, pp. 1-358


Spatial Parallelism of a 3D Finite Difference Velocity-Stress Elastic Wave Propagation Code

Susan E. Minkoff

SIAM J. Sci. Comput. 24, pp. 1-19 (19 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In a three-dimensional isotropic elastic earth, the wave equation solution consists of three velocity components and six stresses. We discretize the partial derivatives using second order in time and fourth order in space staggered finite difference operators. The parallel implementation uses the message passing interface library for platform portability and a spatial decomposition for efficiency. Most of the communication in the code consists of passing subdomain face information to neighboring processors. When the parallel communication is balanced against computation by allocating subdomains of reasonable size, we observe excellent scaled speedup. Allocating subdomains of size $25 \times 25 \times 25$ on each node, we achieve efficiencies of 94\% on 128 processors of an Intel\linebreak Paragon.

GMRES with Deflated Restarting

Ronald B. Morgan

SIAM J. Sci. Comput. 24, pp. 20-37 (18 pages) | Cited 50 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A modification is given of the GMRES iterative method for nonsymmetric systems of linear equations. The new method deflates eigenvalues using Wu and Simon's thick restarting approach [SIAM J. Matrix Anal. Appl., 22 (2000), pp. 602--616]. It has the efficiency of implicit restarting but is simpler and does not have the same numerical concerns. The deflation of small eigenvalues can greatly improve the convergence of restarted GMRES. Also, it is demonstrated that using harmonic Ritz vectors is important because then the whole subspace is a Krylov subspace that contains certain important smaller subspaces.

Optimized Schwarz Methods without Overlap for the Helmholtz Equation

Martin J. Gander, Frédéric Magoulès, and Frédéric Nataf

SIAM J. Sci. Comput. 24, pp. 38-60 (23 pages) | Cited 15 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The classical Schwarz method is a domain decomposition method to solve elliptic partial differential equations in parallel. Convergence is achieved through overlap of the subdomains. We study in this paper a variant of the Schwarz method which converges without overlap for the Helmholtz equation. We show that the key ingredients for such an algorithm are the transmission conditions. We derive optimal transmission conditions which lead to convergence of the algorithm in a finite number of steps. These conditions are, however, nonlocal in nature, and we introduce local approximations which we optimize for performance of the Schwarz method. This leads to an algorithm in the class of optimized Schwarz methods. We present an asymptotic analysis of the optimized Schwarz method for two types of transmission conditions, Robin conditions and transmission conditions with second order tangential derivatives. Numerical results illustrate the effectiveness of the optimized Schwarz method on a model problem and on a problem from industry.

On the Stationary Points of the Squared Distance between Two Ellipses with a Common Focus

G. F. Gronchi

SIAM J. Sci. Comput. 24, pp. 61-80 (20 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we introduce an effective algebraic method for the computation of all the stationary points of the squared distance $d^2$ between a point on one ellipse and a point on a second ellipse with a focus in common with the first one.
This problem comes from celestial mechanics, in which the minima between two elliptic Keplerian orbits are relevant to study the probability of collision; some applications of our algorithm in this field are shown.
This algorithm is based on the use of the fast Fourier transform to obtain the coefficients of the resultant of the two bivariate components of the gradient of $d^2$ with respect to one variable and relies on specific tools in symbolic computation.
An upper bound to the total number of stationary points that we have to expect in this problem is also given; this is done using some tools from algebraic geometry.

A Multiresolution Approach to Regularization of Singular Operators and Fast Summation

Gregory Beylkin and Robert Cramer

SIAM J. Sci. Comput. 24, pp. 81-117 (37 pages) | Cited 14 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Singular and hypersingular operators are ubiquitous in problems of physics, and their use requires a careful numerical interpretation. Although analytical methods for their regularization have long been known, the classical approach does not provide numerical procedures for constructing or applying the regularized operator. We present a multiresolution definition of regularization for integral operators with convolutional kernels which are homogeneous or associated homogeneous functions. We show that our procedure yields the same operator as the classical definition. Moreover, due to the constructive nature of our definition, we provide concise numerical procedures for the construction and application of singular and hypersingular operators. As an application, we present an algorithm for fast computation of discrete sums and briefly discuss several other examples.

A Moving Mesh Method Based on the Geometric Conservation Law

Weiming Cao, Weizhang Huang, and Robert D. Russell

SIAM J. Sci. Comput. 24, pp. 118-142 (25 pages) | Cited 11 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A new adaptive mesh movement strategy is presented, which, unlike many existing moving mesh methods, targets the mesh velocities rather than the mesh coordinates. The mesh velocities are determined in a least squares framework by using the geometric conservation law, specifying a form for the Jacobian determinant of the coordinate transformation defining the mesh, and employing a curl condition. By relating the Jacobian to a monitor function, one is able to directly control the mesh concentration. The geometric conservation law, an identity satisfied by any nonsingular coordinate transformation, is an important tool which has been used for many years in the engineering community to develop cell-volume-preserving finite-volume schemes. It is used here to transform the algebraic expression specifying the Jacobian into an equivalent differential relation which is the key formula for the new mesh movement strategy. It is shown that the resulting method bears a close relation with the Lagrangian method. Advantages of the new approach include the ease of controlling the cell volumes (and therefore mesh adaption) and a theoretical guarantee for existence and nonsingularity of the coordinate transformation. It is shown that the method may suffer from the mesh skewness, a consequence resulting from its close relation with the Lagrangian method. Numerical results are presented to demonstrate various features of the new method.

Accuracy, Resolution, and Stability Properties of a Modified Chebyshev Method

Jodi L. Mead and Rosemary A. Renaut

SIAM J. Sci. Comput. 24, pp. 143-160 (18 pages) | Cited 9 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
While the Chebyshev pseudospectral method provides a spectrally accurate method, integration of partial differential equations with spatial derivatives of order $M$ requires time steps of approximately $O(N^{-2M})$ for stable explicit solvers. Theoretically, time steps may be increased to $O(N^{-M})$ with the use of a parameter, $\alpha$-dependent mapped method introduced by Kosloff and Tal-Ezer [{\em J.\ Comput.\ Phys}., 104 (1993), pp. 457--469]. Our analysis focuses on the utilization of this method for reasonable practical choices for $N$, namely $N \lesssim 30$, as may be needed for two- or three-dimensional modeling. Results presented confirm that spectral accuracy with increasing $N$ is possible both for constant $\alpha$ (Hesthaven, Dinesen, and Lynov [{\em J.\ Comput.\ Phys}., 155 (1999), pp. 287--306]) and for $\alpha$ scaled with $N$, $\alpha$ sufficiently different from $1$ (Don and Solomonoff [{\em SIAM J.\ Sci.\ Comput}., 18 (1997), pp. 1040--1055]). Theoretical bounds, however, show that any realistic choice for $\alpha$, in which both resolution and accuracy considerations are imposed, permits no more than a doubling of the time step for a stable explicit integrator in time, much less than the $O(N)$ improvement claimed by Kosloff and Tal-Ezer. On the other hand, by choosing $\alpha$ carefully, it is possible to improve on the resolution of the Chebyshev method; in particular, one may achieve satisfactory resolution with fewer than $\pi$ points per wavelength. Moreover, this improvement is noted not only for waves with the minimal resolution but also for waves sampled up to about $8$ points per wavelength. Our conclusions are verified by calculation of phase and amplitude errors for numerical solutions of first and second order one-dimensional wave equations. Specifically, while $\alpha$ can be chosen such that the mapped method improves the accuracy and resolution of the Chebyshev method, for practical choices of $N$, it is not possible to achieve both single precision accuracy and gain the advantage of an $O(N^{-M})$ time step.

Fast Convolution for Nonreflecting Boundary Conditions

Christian Lubich and Achim Schädle

SIAM J. Sci. Comput. 24, pp. 161-182 (22 pages) | Cited 37 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Nonreflecting boundary conditions for problems of wave propagation are nonlocal in space and time. While the nonlocality in space can be efficiently handled by Fourier or spherical expansions in special geometries, the arising temporal convolutions still form a computational bottleneck. In the present article, a new algorithm for the evaluation of these convolution integrals is proposed. To compute a temporal convolution over Nt successive time steps, the algorithm requires O(Nt log Nt) operations and O(log Nt) memory. In the numerical examples, this algorithm is used to discretize the Neumann-to-Dirichlet operators arising from the formulation of nonreflecting boundary conditions in rectangular geometries for Schrödinger and wave equations.

Nonlinearly Preconditioned Inexact Newton Algorithms

Xiao-Chuan Cai and David E. Keyes

SIAM J. Sci. Comput. 24, pp. 183-200 (18 pages) | Cited 22 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Inexact Newton algorithms are commonly used for solving large sparse nonlinear system of equations $F(u^{\ast})=0$ arising, for example, from the discretization of partial differential equations. Even with global strategies such as linesearch or trust region, the methods often stagnate at local minima of $\|F\|$, especially for problems with unbalanced nonlinearities, because the methods do not have built-in machinery to deal with the unbalanced nonlinearities. To find the same solution $u^{\ast}$, one may want to solve instead an equivalent nonlinearly preconditioned system ${\cal F}(u^{\ast})=0$ whose nonlinearities are more balanced. In this paper, we propose and study a nonlinear additive Schwarz-based parallel nonlinear preconditioner and show numerically that the new method converges well even for some difficult problems, such as high Reynolds number flows, where a traditional inexact Newton method fails.

Adjusting the Rayleigh Quotient in Semiorthogonal Lanczos Methods

G. W. Stewart

SIAM J. Sci. Comput. 24, pp. 201-207 (7 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In a semiorthogonal Lanczos algorithm, the orthogonality of the Lanczos vectors is allowed to deteriorate to roughly the square root of the rounding unit, after which the current vectors are reorthogonalized. A theorem of Simon [Linear Algebra Appl., 61 (1984), pp. 101--132] shows that the Rayleigh quotient---i.e., the tridiagonal matrix produced by the Lanczos recursion---contains fully accurate approximations to the Ritz values in spite of the lack of orthogonality. Unfortunately, the same lack of orthogonality can cause the Ritz vectors to fail to converge. It also makes the classical estimate for the residual norm misleadingly small. In this paper we show how to adjust the Rayleigh quotient to overcome this problem.

Interface Tracking for Axisymmetric Flows

James Glimm, John W. Grove, and Yongmin Zhang

SIAM J. Sci. Comput. 24, pp. 208-236 (29 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A front tracking method for inviscid gas dynamics is presented. The key constructions and algorithms used in the code are described and the interrelations between shock capturing, interface dynamics, computational geometry, grid construction, and parallelism are discussed for the code as a whole. Validation is carried out by comparing the single mode bubble velocity for Rayleigh--Taylor instability with theoretical models and experimental results. The calculations are validated by mesh refinement studies and by the comparison of the asymptotic limit of the minimum radius $r_{\rm min} \rightarrow \infty$ to a pure planar computation in two dimensions.

A Preconditioner for the Steady-State Navier--Stokes Equations

David Kay, Daniel Loghin, and Andrew Wathen

SIAM J. Sci. Comput. 24, pp. 237-256 (20 pages) | Cited 57 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We present a new method for solving the sparse linear system of equations arising from the discretization of the linearized steady-state Navier--Stokes equations (also known as the Oseen equations). The solver is an iterative method of Krylov subspace type for which we devise a preconditioner through a heuristic argument based on the fundamental solution tensor for the Oseen operator. The preconditioner may also be conceived through a weaker heuristic argument involving differential operators. Computations indicate that convergence for the preconditioned discrete Oseen problem is only mildly dependent on the viscosity (inverse Reynolds number) and, most importantly, that the number of iterations does not grow as the mesh size is reduced. Indeed, since the preconditioner is motivated through analysis of continuous operators, the number of iterations decreases for smaller mesh size which accords with better approximation of these operators.

A New Error Estimate of the Fast Gauss Transform

B. J. C. Baxter and George Roussos

SIAM J. Sci. Comput. 24, pp. 257-259 (3 pages) | Cited 9 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The fast Gauss transform of L. Greengard and J. Strain [SIAM J. Sci. Statist. Comput., 12 (1991), pp. 79--94] reduces the computational complexity of the evaluation of the sum of N Gaussians at M points in d-dimensional space from ${\cal O}(MN)$ to ${\cal O}(M+N)$ floating-point operations. In this note, we provide numerical evidence that the error estimate of Lemma 2.1 in [SIAM J. Sci. Statist. Comput., 12 (1991), pp. 79--94] is erroneous and then proceed to calculate a replacement error estimate for the fast Gauss transform, incorporating an improved upper bound for Hermite functions.

On Two Variants of an Algebraic Wavelet Preconditioner

Tony F. Chan and Ke Chen

SIAM J. Sci. Comput. 24, pp. 260-283 (24 pages) | Cited 8 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A recursive method of constructing preconditioning matrices for the nonsymmetric stiffness matrix in a wavelet basis is proposed for solving a class of integral and differential equations. It is based on a level-by-level application of the wavelet scales decoupling the different wavelet levels in a matrix form just as in the well-known nonstandard form. The result is a powerful iterative method with built-in preconditioning leading to two specific algebraic multilevel iteration algorithms: one with an exact Schur preconditioning and the other with an approximate Schur preconditioning. Numerical examples are presented to illustrate the efficiency of the new algorithms.

A Parallel Implementation of the Nonsymmetric QR Algorithm for Distributed Memory Architectures

Greg Henry, David Watkins, and Jack Dongarra

SIAM J. Sci. Comput. 24, pp. 284-311 (28 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
One approach to solving the nonsymmetric eigenvalue problem in parallel is to parallelize the QR algorithm. Not long ago, this was widely considered to be a hopeless task. Recent efforts have led to significant advances, although the methods proposed up to now have suffered from scalability problems. This paper discusses an approach to parallelizing the QR algorithm that greatly improves scalability. A theoretical analysis indicates that the algorithm is ultimately not scalable, but the nonscalability does not become evident until the matrix dimension is enormous. Experiments on the Intel Paragon system, the IBM SP2 supercomputer, the SGI Origin 2000, and the Intel ASCI Option Red supercomputer are reported.

An Inverse Free Preconditioned Krylov Subspace Method for Symmetric Generalized Eigenvalue Problems

Gene H. Golub and Qiang Ye

SIAM J. Sci. Comput. 24, pp. 312-334 (23 pages) | Cited 9 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we present an inverse free Krylov subspace method for finding some extreme eigenvalues of the symmetric definite generalized eigenvalue problem $Ax = \lambda B x$. The basic method takes a form of inner-outer iterations and involves no inversion of B or any shift-and-invert matrix $A-\lambda_0 B$. A convergence analysis is presented that leads to a preconditioning scheme for accelerating convergence through some equivalent transformations of the eigenvalue problem. Numerical examples are given to illustrate the convergence properties and to demonstrate the competitiveness of the method.

Stabilized Explicit-Implicit Domain Decomposition Methods for the Numerical Solution of Parabolic Equations

Yu Zhuang and Xian-He Sun

SIAM J. Sci. Comput. 24, pp. 335-358 (24 pages) | Cited 12 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We report a class of stabilized explicit-implicit domain decomposition (SEIDD) methods for the numerical solution of parabolic equations. Explicit-implicit domain decomposition (EIDD) methods are globally noniterative, nonoverlapping domain decomposition methods, which, when compared with Schwarz-algorithm-based parabolic solvers, are computationally and communicationally efficient for each simulation time step but suffer from small time step size restrictions. By adding a stabilization step to EIDD, the SEIDD methods retain the time-stepwise efficiency in computation and communication of the EIDD methods but exhibit much better numerical stability. Three SEIDD algorithms are presented in this paper, which are experimentally tested to show excellent stability, computation and communication efficiencies, and high parallel speedup and scalability.
Close

Your Recent History Recent History Help

Recently Viewed

  1. Iterative Solution of Cyclically Reduced Systems Arising from Discretization of the Three-Dimensiona...
  2. Nonoscillatory Central Schemes for Multidimensional Hyperbolic Conservation Laws
  3. A Fast Multilevel Algorithm for the Solution of Nonlinear Systems of Conductive-Radiative Heat Trans...
  4. Parallel Newton--Krylov--Schwarz Algorithms for the Transonic Full Potential Equation
  5. Numerical Solution of the Riemann Problem for Two-Dimensional Gas Dynamics
© SIAM By using SIAM Publications Online you agree to abide by the Terms and Conditions of Use. Banner art adapted from a figure by Hinke M. Osinga and Bernd Krauskopf (University of Bristol, UK).
close