SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

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

Search Issue | RSS Feeds RSS
Previous Issue

2001

Volume 22, Issue 6, pp. 1905-2281


Structure-Preserving Methods for Computing Eigenpairs of Large Sparse Skew-Hamiltonian/Hamiltonian Pencils

Volker Mehrmann and David Watkins

SIAM J. Sci. Comput. 22, pp. 1905-1925 (21 pages) | Cited 27 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We study large, sparse generalized eigenvalue problems for matrix pencils, where one of the matrices is Hamiltonian and the other is skew-Hamiltonian. Problems of this form arise in the numerical simulation of elastic deformation of anisotropic materials, in structural mechanics, and in the linear quadratic control problem for partial differential equations. We develop a structure-preserving skew-Hamiltonian, isotropic, implicitly restarted shift-and-invert Arnoldi algorithm (SHIRA). Several numerical examples demonstrate the superiority of SHIRA over a competing unstructured method.

Comparison of the Nodal Integral Method and Nonstandard Finite-Difference Schemes for the Fisher Equation

Rizwan Uddin

SIAM J. Sci. Comput. 22, pp. 1926-1942 (17 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The relationship between the so-called nonstandard finite-difference schemes and the nodal integral method (NIM) is investigated using the Fisher equation as a model problem. Exact and best finite-difference schemes are reviewed first. Next, the NIM for the Fisher equation is developed. It is shown that the NIM leads to a nonstandard evaluation of the derivatives. Moreover, the resulting scheme possesses the desirable characteristics of the nonstandard finite difference schemes, such as the nonlocal evaluation of the nonlinear terms. Thus, the NIM provides a systematic framework to obtain schemes similar to the best finite-difference schemes. Numerical results for a propagating front problem show that the NIM can capture the shape and speed of the front very accurately. Results also show that the best finite-difference scheme is stable for large grid sizes but only at the cost of inaccuracy in the front propagation speed. Additional results are obtained using the NIM for symmetric and asymmetric initial conditions. These describe the interaction of two fronts of advantageous genes that approach each other and merge. It is noted that the traveling fronts evolving from certain asymmetric initial conditions become indistinguishable from those that evolve from symmetric initial conditions.

Fast Finite Volume Simulation of 3D Electromagnetic Problems with Highly Discontinuous Coefficients

E. Haber and U. M. Ascher

SIAM J. Sci. Comput. 22, pp. 1943-1961 (19 pages) | Cited 32 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We consider solving three-dimensional electromagnetic problems in parameter regimes where the quasi-static approximation applies and the permeability, permittivity, and conductivity may vary significantly. The difficulties encountered include handling solution discontinuities across interfaces and accelerating convergence of traditional iterative methods for the solution of the linear systems of algebraic equations that arise when discretizing Maxwell's equations in the frequency domain.
The present article extends methods we proposed earlier for constant permeability [E. Haber, U. Ascher, D. Aruliah, and D. Oldenburg, J. Comput. Phys., 163 (2000), pp. 150--171; D. Aruliah, U. Ascher, E. Haber, and D. Oldenburg, Math. Models Methods Appl. Sci., to appear.] to handle also problems in which the permeability is variable and may contain significant jump discontinuities. In order to address the problem of slow convergence we reformulate Maxwell's equations in terms of potentials, applying a Helmholtz decomposition to either the electric field or the magnetic field. The null space of the curl operators can then be annihilated by adding a stabilizing term, using a gauge condition, and thus obtaining a strongly elliptic differential operator. A staggered grid finite volume discretization is subsequently applied to the reformulated PDE system. This scheme works well for sources of various types, even in the presence of strong material discontinuities in both conductivity and permeability. The resulting discrete system is amenable to fast convergence of ILU-preconditioned Krylov methods.
We test our method using several numerical examples and demonstrate its robust efficiency. We also compare it to the classical Yee method using similar iterative techniques for the resulting algebraic system, and we show that our method is significantly faster, especially for electric sources.

Exact Prediction of QR Fill-In by Row-Merge Trees

Suely Oliveira

SIAM J. Sci. Comput. 22, pp. 1962-1973 (12 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Row-merge trees for forming the QR factorization of a sparse matrix A are closely related to elimination trees for the Cholesky factorization of AT A. Row-merge trees predict the exact fill-in (assuming no numerical cancellation) provided A satisfies the strong Hall property, but overestimate the fill-in in general. However, here a fast and simple postprocessing step for row-merge trees is presented that predicts the exact fill-in for sparse QR factorization using Householder reflectors for general matrices.

On Schwarz Alternating Methods for the Incompressible Navier--Stokes Equations

S. H. Lui

SIAM J. Sci. Comput. 22, pp. 1974-1986 (13 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The Schwarz alternating method can be used to solve linear elliptic boundary value problems on domains which consist of two or more overlapping subdomains. The solution is approximated by an infinite sequence of functions which result from solving a sequence of elliptic boundary value problems in each of the subdomains.
This paper considers four Schwarz alternating methods for the N-dimensional, steady, viscous, incompressible Navier--Stokes equations, $N \leq 4$. It is shown that the Schwarz sequences converge to the true solution provided that the Reynolds number is sufficiently small.

Algebraic Two-Level Preconditioners for the Schur Complement Method

L. M. Carvalho, L. Giraud, and P. Le Tallec

SIAM J. Sci. Comput. 22, pp. 1987-2005 (19 pages) | Cited 12 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The solution of elliptic problems is challenging on parallel distributed memory computers since their Green's functions are global. To address this issue, we present a set of preconditioners for the Schur complement domain decomposition method. They implement a global coupling mechanism, through coarse-space components, similar to the one proposed in [Bramble, Pasciak, and Shatz, Math. Comp., 47 (1986), pp. 103--134]. The definition of the coarse-space components is algebraic; they are defined using the mesh partitioning information and simple interpolation operators. These preconditioners are implemented on distributed memory computers without introducing any new global synchronization in the preconditioned conjugate gradient iteration. The numerical and parallel scalability of those preconditioners are illustrated on two-dimensional model examples that have anisotropy and/or discontinuity phenomena.

Asymptotic Analysis of the Laminar Viscous Flow Over a Porous Bed

Willi Jäger, Andro Mikelic, and Nicolas Neuss

SIAM J. Sci. Comput. 22, pp. 2006-2028 (23 pages) | Cited 17 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We consider the laminar viscous channel flow over a porous surface. The size of the pores is much smaller than the size of the channel, and it is important to determine the effective boundary conditions at the porous surface. We study the corresponding boundary layers, and, by a rigorous asymptotic expansion, we obtain Saffman's modification of the interface condition observed by Beavers and Joseph. The effective coefficient in the law is determined through an auxiliary boundary-layer type problem, whose computational and modeling aspects are discussed in detail. Furthermore, the approximation errors for the velocity and for the effective mass flow are given as powers of the characteristic pore size $\ep$. Finally, we give the interface condition linking the effective pressure fields in the porous medium and in the channel, and we determine the jump of the effective pressures explicitly.

Special Bilinear Quadrilateral Elements For Locally Refined Finite Element Grids

Weigang Wang

SIAM J. Sci. Comput. 22, pp. 2029-2050 (22 pages) | Cited 3 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A set of special bilinear quadrilateral elements are developed to handle irregular nodes on the interface between different refinement levels in 1-irregular meshes for the adaptive finite element method. With this alternative approach, irregular nodes are eliminated, and at the same time the C0 continuity at interelements is strictly ensured. The new elements simplify the refinement procedure and circumvent many inconveniences in local refinement with linear quadrilateral elements. Solutions to problems of a heat conduction, a reaction-diffusion, and incompressible and viscous flows past a circular cylinder are respectively presented to test the effectiveness of special elements. They demonstrate good accuracy and a potential applicability to more complex problems.

Constraint Partitioning for Stability in Path-Constrained Dynamic Optimization Problems

Soumyendu Raha and Linda R. Petzold

SIAM J. Sci. Comput. 22, pp. 2051-2074 (24 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper an algorithm for extracting a stable differential-algebraic subsystem from a path-constrained dynamical system is proposed. The subsystem may be integrated directly by a differential-algebraic system integrator to evaluate constraints in shooting- or multiple shooting-type direct methods for solving path-constrained dynamic optimization problems. The algorithm appends algebraic constraints to the unconstrained ordinary differential equation subsystem based on a stability estimate for the resulting differential-algebraic system. The logarithmic norm is used to compute a stability estimate for index 1 and index 2 subsystems. The working of the algorithm is illustrated with examples.

Two Element-by-Element Iterative Solutions for Shallow Water Equations

C. C. Fang and Tony W. H. Sheu

SIAM J. Sci. Comput. 22, pp. 2075-2092 (18 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we apply the generalized Taylor--Galerkin finite element model to simulate bore wave propagation in a domain of two dimensions. For stability and accuracy reasons, we generalize the model through the introduction of four free parameters. One set of parameters is rigorously determined to obtain the high-order finite element solution. The other set of free parameters is determined from the underlying discrete maximum principle to obtain the monotonic solutions. The resulting two models are used in combination through the flux correct transport technique of Zalesak, thereby constructing a finite element model which has the ability to capture hydraulic discontinuities. In addition, this paper highlights the implementation of two Krylov subspace iterative solvers, namely, the bi-conjugate gradient stabilized (Bi-CGSTAB) and the generalized minimum residual (GMRES) methods. For the sake of comparison, the multifrontal direct solver is also considered. The performance characteristics of the investigated solvers are assessed using results of a standard test widely used as a benchmark in hydraulic modeling. Based on numerical results, it is shown that the present finite element method can render the technique suitable for solving shallow water equations with sharply varying solution profiles. Also, the GMRES solver is shown to have a much better convergence rate than the Bi-CGSTAB solver, thereby saving much computing time compared to the multifrontal solver.

Coulomb Interactions on Planar Structures: Inverting the Square Root of the Laplacian

Zydrunas Gimbutas, Leslie Greengard, and Michael Minion

SIAM J. Sci. Comput. 22, pp. 2093-2108 (16 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We present an adaptive fast multipole method for inverting the square root of the Laplacian in two dimensions. Solving this problem is the dominant computational cost in many applications arising in electrical engineering, geophysical fluid dynamics, and the study of thin films. It corresponds to the evaluation of the field induced by a planar distribution of charge or vorticity. Our algorithm is direct and assumes only that the source distribution is discretized using an adaptive quad-tree. The amount of work grows linearly with the number of mesh points.

CGNR Is an Error Reducing Algorithm

Chunguang Li

SIAM J. Sci. Comput. 22, pp. 2109-2112 (4 pages)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The CGNR algorithm is a robust algorithm solving large nonsymmetric linear systems. In this short note, we show that the norm of the error vector of CGNR approximation decreases monotonically during its iteration. Few iterative algorithms have this error reducing property. The result indicates that the simple CGNR algorithm is still a potential one in practical use.

Fast Bit-Reversals on Uniprocessors and Shared-Memory Multiprocessors

Zhao Zhang and Xiaodong Zhang

SIAM J. Sci. Comput. 22, pp. 2113-2134 (22 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we examine different methods using techniques of blocking, buffering, and padding for efficient implementations of bit-reversals. We evaluate the merits and limits of each technique and its application and architecture-dependent conditions for developing cache-optimal methods. Besides testing the methods on different uniprocessors, we conducted both simulation and measurements on two commercial symmetric multiprocessors (SMP) to provide architectural insights into the methods and their implementations. We present two contributions in this paper: (1) Our integrated blocking methods, which match cache associativity and translation-lookaside buffer (TLB) cache size and which fully use the available registers, are cache-optimal and fast. (2) We show that our padding methods outperform other software-oriented methods, and we believe they are the fastest in terms of minimizing both CPU and memory access cycles. Since the padding methods are almost independent of hardware, they could be widely used on many uniprocessor workstations and multiprocessors.

Checkpointing Schemes for Adjoint Codes: Application to the Meteorological Model Meso-NH

I. Charpentier

SIAM J. Sci. Comput. 22, pp. 2135-2151 (17 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The adjoint code of a nonlinear computer model calculates gradients along a trajectory that has to be known at integration time. When the storage of the whole trajectory requires too large an amount of memory, the calculation of the adjoint code is split and is done part by part from restart points called checkpoints. Griewank proposed a checkpointing method named Revolve, which provides an optimal logarithmic behavior with respect to time and memory requirement. In this work, some checkpointing schedules are proposed. Some of them correspond to special cases of Revolve.
The user's preference is essential to choose between time and memory requirements. This is a key point for adjoint codes of temporal models such as the meteorological model Meso-NH that may be used for weather forecasts. When the computational time is the top priority, a particular checkpointing scheme allows computation of the adjoint code with at most one extra integration of the model. The memory requirement behaves then as the square root of the number of iterations of the model. Checkpointing schemes are tested on adjoint simulations of Meso-NH.

A Semi-Lagrangian Finite Volume Method for Newtonian Contraction Flows

T. N. Phillips and A. J. Williams

SIAM J. Sci. Comput. 22, pp. 2152-2177 (26 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A new finite volume method for solving the incompressible Navier--Stokes equations is presented. The main features of this method are the location of the velocity components and pressure on different staggered grids and a semi-Lagrangian method for the treatment of convection. An interpolation procedure based on area-weighting is used for the convection part of the computation. The method is applied to flow through a constricted channel, and results are obtained for Reynolds numbers, based on half the flow rate, up to 1000. The behavior of the vortex in the salient corner is investigated qualitatively and quantitatively, and excellent agreement is found with the numerical results of Dennis and Smith [ Proc. Roy. Soc. London A, 372 (1980), pp. 393--414] and the asymptotic theory of Smith [J. Fluid Mech., 90 (1979), pp. 725--754].

An $\cal O(N)$ Level Set Method for Eikonal Equations

Seongjai Kim

SIAM J. Sci. Comput. 22, pp. 2178-2193 (16 pages) | Cited 15 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A propagating interface can develop corners and discontinuities as it advances. Level set algorithms have been extensively applied for the problems in which the solution has advancing fronts. One of the most popular level set algorithms is the so-called {fast marching method} (FMM), which requires total $\cal O(N\log_2N)$ operations, where N is the number of grid points. The article is concerned with the development of an $\cal O(N)$ level set algorithm called the group marching method (GMM). The new method is based on the narrow band approach as in the FMM. However, it is incorporating a correction-by-iteration strategy to advance a group of grid points at a time, rather than sorting the solution in the narrow band to march forward a single grid point. After selecting a group of grid points appropriately, the GMM advances the group in two iterations for the cost of slightly larger than one iteration. Numerical results are presented to show the efficiency of the method, applied to the eikonal equation in two and three dimensions.

A Scalable Parallel Algorithm for Incomplete Factor Preconditioning

David Hysom and Alex Pothen

SIAM J. Sci. Comput. 22, pp. 2194-2215 (22 pages) | Cited 24 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We describe a parallel algorithm for computing incomplete factor (ILU) preconditioners. The algorithm attains a high degree of parallelism through graph partitioning and a two-level ordering strategy. Both the subdomains and the nodes within each subdomain are ordered to preserve concurrency. We show through an algorithmic analysis and through computational results that this algorithm is scalable. Experimental results include timings on three parallel platforms for problems with up to 20 million unknowns running on up to 216 processors. The resulting preconditioned Krylov solvers have the desirable property that the number of iterations required for convergence is insensitive to the number of processors.

Accurate Computations on Inertial Manifolds

M. S. Jolly, R. Rosa, and R. Temam

SIAM J. Sci. Comput. 22, pp. 2216-2238 (23 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
An algorithm is implemented for the computation of inertial manifolds to arbitrary accuracy. It is applied to the Kuramoto--Sivashinsky equation to recover the high wave number components of elements on the global attractor. The computational results demonstrate that this algorithm is significantly more accurate than previous methods, when compared after just a few iterations. An essential feature of the algorithm is that it does not require that the entire manifold be constructed. Hence the algorithm is particularly suited for computing trajectories on the manifold, and its computational complexity is independent of the dimension.

Initial Values for the Inverse Toeplitz Eigenvalue Problem

Dirk Laurie

SIAM J. Sci. Comput. 22, pp. 2239-2255 (17 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We show how to find a symmetric Toeplitz matrix whose eigenvalues are close to a specified target set of n real numbers.

High Resolution Schemes for Conservation Laws with Locally Varying Time Steps

Clint Dawson and Robert Kirby

SIAM J. Sci. Comput. 22, pp. 2256-2281 (26 pages) | Cited 22 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We develop upwind methods which use limited high resolution corrections in the spatial discretization and local time stepping for forward Euler and second order time discretizations. $L^\infty$ stability is proven for both time stepping schemes for problems in one space dimension. These methods are restricted by a local CFL condition rather than the traditional global CFL condition, allowing local time refinement to be coupled with local spatial refinement. Numerical evidence demonstrates the stability and accuracy of the methods for problems in both one and two space dimensions.
Close

close