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

2004

Volume 25, Issue 6, pp. 1837-2187


Encapsulating Multiple Communication-Cost Metrics in Partitioning Sparse Rectangular Matrices for Parallel Matrix-Vector Multiplies

Bora Uçar and Cevdet Aykanat

SIAM J. Sci. Comput. 25, pp. 1837-1859 (23 pages) | Cited 10 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This paper addresses the problem of one-dimensional partitioning of structurally unsymmetric square and rectangular sparse matrices for parallel matrix-vector and matrix-transpose-vector multiplies. The objectiveis to minimize the communication cost while maintaining the balance on computational loads of processors. Most of the existing partitioning models consider only the total message volume hoping that minimizing this communication-cost metric is likely to reduce other metrics. However, the total message latency (start-up time) may be more important than the total message volume. Furthermore, the maximum message volume and latency handled by a single processor are also important metrics. We propose a two-phase approach that encapsulates all these four communication-cost metrics. The objective in the first phase is to minimize the total message volume while maintaining the computational-load balance. The objective in the second phase is to encapsulate the remaining three communication-cost metrics. We propose communication-hypergraph and partitioning models for the second phase. We then present several methods for partitioning communication hypergraphs. Experiments on a wide range of test matrices show that the proposed approach yields very effective partitioning results. A parallel implementation on a PC cluster verifies that the theoretical improvements shown by partitioning results hold in practice.

Permuting Sparse Rectangular Matrices into Block-Diagonal Form

Cevdet Aykanat, Ali Pinar, and Ümit V. Çatalyürek

SIAM J. Sci. Comput. 25, pp. 1860-1879 (20 pages) | Cited 23 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We investigate the problem of permuting a sparse rectangular matrix into block-diagonal form. Block-diagonal form of a matrix grants an inherent parallelism for solving the deriving problem, as recently investigated in the context of mathematical programming, LU factorization, and QR factorization. To represent the nonzero structure of a matrix, we propose bipartite graph and hypergraph models that reduce the permutation problem to those of graph partitioning by vertex separator and hypergraph partitioning, respectively. Our experiments on a wide range of matrices, using the state-of-the-art graph and hypergraph partitioning tools MeTiS and PaToH\@, revealed that the proposed methods yield very effective solutions both in terms of solution quality and runtime.

Dynamic Database Generation for Efficient Calculation of Stellarator Plasma Equilibria

S. P. Hirshman, L. A. Berry, and S. Jesse

SIAM J. Sci. Comput. 25, pp. 1880-1895 (16 pages)

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
A numerical scheme has been developed and tested that significantly reduces (by a factor of 30, typically) both the computation time and storage requirements for external magnetic field calculations embedded in a three-dimensional iterative magnetohydrodynamic equilibrium calculation. The equilibrium is computed from external magnetic fields that are calculated numerically using the Biot--Savart law applied to a complex set of coils. These fields must be evaluated on the plasma boundary, which evolves as the equilibrium converges. The improved efficiency of the present method results from dynamically building and storing a database of magnetic field components on a grid that is determined by the changing plasma boundary itself. When possible, field values on the plasma boundary for a particular iteration are interpolated from these grid values, thus eliminating the need for further time-consuming Biot--Savart calculations. New grid points are added to the database only when the boundary moves outside the spatial range of previously computed interpolation values. This dynamical method is particularly efficient when the equilibrium calculation is embedded in an optimization loop, where evolution of the coil set requires many reevaluations of the database.

Adaptive Smoothed Aggregation ($\alpha$SA)

M. Brezina, R. Falgout, S. MacLachlan, T. Manteuffel, S. McCormick, and J. Ruge

SIAM J. Sci. Comput. 25, pp. 1896-1920 (25 pages) | Cited 34 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Substantial effort has been focused over the last two decades on developing multilevel iterative methods capable of solving the large linear systems encountered in engineering practice. These systems often arise from discretizing partial differential equations over unstructured meshes, and the particular parameters or geometry of the physical problem being discretized may be unavailable to the solver. Algebraic multigrid (AMG) and multilevel domain decomposition methods of algebraic type have been of particular interest in this context because of their promises of optimal performance without the need for explicit knowledge of the problem geometry. These methods construct a hierarchy of coarse problems based on the linear system itself and on certain assumptions about the smooth components of the error. For smoothed aggregation (SA) methods applied to discretizations of elliptic problems, these assumptions typically consist of knowledge of the near-nullspace of the weak form. This paper introduces an extension of the SA method in which good convergence properties are achieved in situations where explicit knowledge of the near-nullspace components is unavailable. This extension is accomplished by using the method itself to determine near-nullspace components and adjusting the coarsening processes accordingly.

A Broyden Rank p+1 Update Continuation Method with Subspace Iteration

T. L. van Noorden, S. M. Verduyn Lunel, and A. Bliek

SIAM J. Sci. Comput. 25, pp. 1921-1940 (20 pages) | Cited 1 time

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper we present an efficient branch-following procedure that can be used not only to compute branches of periodic solutions of periodically forced dynamical systems but also to determine the stability of the periodic solutions. The procedure combines Broyden's method with a subspace iteration method to determine the dominant eigenvalues. The method has connections with the hybrid Newton--Picard methods developed by Lust et al. in [SIAM J. Sci. Comput., 19 (1998) pp. 1188--1209]. A convergence analysis of the procedure is presented. The method is applied to the computation of periodic states of a reverse flow reactor, and its performance is compared with two variants of the hybrid Newton--Picard method.

Discontinuous Galerkin BGK Method for Viscous Flow Equations: One-Dimensional Systems

Kun Xu

SIAM J. Sci. Comput. 25, pp. 1941-1963 (23 pages) | Cited 6 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This paper is about the construction of a BGK Navier--Stokes (BGK-NS) solver in the discontinuous Galerkin (DG) framework. Since in the DG formulation the conservative variables and their slopes can be updated simultaneously, the flow evolution in each element involves only the flow variables in the nearest neighboring cells. Instead of using the semidiscrete approach in the Runge--Kutta discontinuous Galerkin (RKDG) method, the current DG-BGK method integrates the governing equations in time as well. Due to the coupling of advection and dissipative terms in the gas-kinetic formulation, the DG-BGK method solves the viscous governing equations directly. Numerical examples for the one-dimensional compressible Navier--Stokes solutions will be presented.

Robust Variance Reduction for Random Walk Methods

Gang Zou and Robert D. Skeel

SIAM J. Sci. Comput. 25, pp. 1964-1981 (18 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Random walk methods are effective for solving linear partial differential equations in many dimensions, especially those involving complex geometries. They are based on an equivalence given by a Feynman--Kac formula between an expectation of a functional of a stochastic process and the solution at a point of a partial differential equation. The drawback is that the error is proportional only to the square root of the reciprocal of the number of trials. Efficiency depends critically on variance reduction. A general strategy for doing this in the case of stochastic differential equations is proposed by Milstein. The idea is to introduce a bias in the drift term and to exactly compensate for this by unequal weighting of the trials. There is an optimal bias defined in terms of the solution of the partial differential equation which reduces the variance to zero. In practice, an approximation is used. This idea has been tested under the name "biased Brownian dynamics" on the problem of calculating rate constants for diffusion-limited reactions. The approach is successful in some cases but is less successful in more difficult cases due to the occasional occurrence of a well-above-average weight. Proposed and tested here is a weight control algorithm, which greatly enhances the effectiveness of biased Brownian dynamics.

A Boundary Condition--Capturing Multigrid Approach to Irregular Boundary Problems

Justin W. L. Wan and Xu-Dong Liu

SIAM J. Sci. Comput. 25, pp. 1982-2003 (22 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We propose a geometric multigrid method for solving linear systems arising from irregular boundary problems involving multiple interfaces in two and three dimensions. In this method, we adopt a matrix-free approach; i.e., we do not form the fine grid matrix explicitly and we never form nor store the coarse grid matrices, as many other robust multigrid methods do. The main idea is to construct an accurate interpolation which captures the correct boundary conditions at the interfaces via a level set function. Numerical results are given to compare our multigrid method with black box and algebraic multigrid methods.

Numerical Simulation of the Smoluchowski Coagulation Equation

Francis Filbet and Philippe Laurençot

SIAM J. Sci. Comput. 25, pp. 2004-2028 (25 pages) | Cited 20 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we develop a numerical scheme for the Smoluchowski coagulation equation, which relies on a conservative formulation and a finite volume approach. Several numerical simulations are performed to test the validity of the scheme and the expected behavior of the model. In particular the gelation phenomenon and the long time behavior of the solution are numerically studied.

Analysis of Preconditioners for Saddle-Point Problems

D. Loghin and A. J. Wathen

SIAM J. Sci. Comput. 25, pp. 2029-2049 (21 pages) | Cited 12 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Mixed finite element formulations give rise to large, sparse, block linear systems of equations, the solution of which is often sought via a preconditioned iterative technique. In this work we present a general analysis of block-preconditioners based on the stability conditions inherited from the formulation of the finite element method (the Babuska--Brezzi, or inf-sup, conditions). The analysis is motivated by the notions of norm-equivalence and field-of-values-equivalence of matrices. In particular, we give sufficient conditions for diagonal and triangular block-preconditioners to be norm- and field-of-values-equivalent to the system matrix.

A Fast and Stable Well-Balanced Scheme with Hydrostatic Reconstruction for Shallow Water Flows

Emmanuel Audusse, François Bouchut, Marie-Odile Bristeau, Rupert Klein, and Benoît Perthame

SIAM J. Sci. Comput. 25, pp. 2050-2065 (16 pages) | Cited 115 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We consider the Saint-Venant system for shallow water flows, with nonflat bottom. It is a hyperbolic system of conservation laws that approximately describes various geophysical flows, such as rivers, coastal areas, and oceans when completed with a Coriolis term, or granular flows when completed with friction. Numerical approximate solutions to this system may be generated using conservative finite volume methods, which are known to properly handle shocks and contact discontinuities. However, in general these schemes are known to be quite inaccurate for near steady states, as the structure of their numerical truncation errors is generally not compatible with exact physical steady state conditions. This difficulty can be overcome by using the so-called well-balanced schemes. We describe a general strategy, based on a local hydrostatic reconstruction, that allows us to derive a well-balanced scheme from any given numerical flux for the homogeneous problem. Whenever the initial solver satisfies some classical stability properties, it yields a simple and fast well-balanced scheme that preserves the nonnegativity of the water height and satisfies a semidiscrete entropy inequality.

Simulations of the Whirling Instability by the Immersed Boundary Method

Sookkyung Lim and Charles S. Peskin

SIAM J. Sci. Comput. 25, pp. 2066-2083 (18 pages) | Cited 23 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
When an elastic filament spins in a viscous incompressible fluid it may undergo a whirling instability, as studied asymptotically by Wolgemuth, Powers, and Goldstein [Phys. Rev. Lett., 84 (2000), pp. 16--23]. We use the immersed boundary (IB) method to study the interaction between the elastic filament and the surrounding viscous fluid as governed by the incompressible Navier--Stokes equations. This allows the study of the whirling motion when the shape of the filament is very different from the unperturbed straight state.

An Automated Multilevel Substructuring Method for Eigenspace Computation in Linear Elastodynamics

Jeffrey K. Bennighof and R. B. Lehoucq

SIAM J. Sci. Comput. 25, pp. 2084-2106 (23 pages) | Cited 20 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
We present an automated multilevel substructuring (AMLS) method for eigenvalue computations in linear elastodynamics in a variational and algebraic setting. AMLS first recursively partitions the domain of the PDE into a hierarchy of subdomains. Then AMLS recursively generates a subspace for approximating the eigenvectors associated with the smallest eigenvalues by computing partial eigensolutions associated with the subdomains and the interfaces between them. We remark that although we present AMLS for linear elastodynamics, our formulation is abstract and applies to generic H1-elliptic bilinear forms.
In the variational formulation, we define an interface mass operator that is consistent with the treatment of elastic properties by the familiar Steklov--Poincaré operator. With this interface mass operator, all of the subdomain and interface eigenvalue problems in AMLS become orthogonal projections of the global eigenvalue problem onto a hierarchy of subspaces. Convergence of AMLS is determined in the continuous setting by the truncation of these eigenspaces, independent of other discretization schemes.
The goal of AMLS, in the algebraic setting, is to achieve a high level of dimensional reduction, locally and inexpensively, while balancing the errors associated with truncation and the finite element discretization. This is accomplished by matching the mesh-independent AMLS truncation error with the finite element discretization error. Our report ends with numerical experiments that demonstrate the effectiveness of AMLS on a model problem and an industrial problem.

Linear versus Nonlinear Dimensionality Reduction of High-Dimensional Dynamical Systems

Nejib Smaoui

SIAM J. Sci. Comput. 25, pp. 2107-2125 (19 pages) | Cited 2 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
Two techniques for dimensionality reduction of high-dimensional dynamical systems are presented. The first is based on Karhunen--Loève (K-L) analysis and the second on autoassociative neural networks (ANNs). First, we analyze the dynamics of two partial differential equations, namely, the one-dimensional (1-d) Kuramoto--Sivashinsky (K-S) equation and the two-dimensional (2-d) Navier--Stokes (N-S) equations. For the 1-d K-S equation, one particular dynamical behavior, represented by a heteroclinic connection in phase space, is investigated. As for the 2-d N-S equations, a quasi-periodic behavior is examined. Coherent structures of both dynamics were extracted spanning linear subspaces with minimum information loss. Then we obtain systems of ordinary differential equations based on sophisticated (K-L) Galerkin-type approximation capturing the dynamics of the attractors of both regimes residing on linear manifolds. Using the K-L data coefficients as inputs to autoassociative neural networks, we are able to reproduce the dynamics of the attractors which now reside on nonlinear manifolds with fewer dimensions than those previously obtained.

A Minimal Residual Interpolation Method for Linear Equations with Multiple Right-Hand Sides

Per Lötstedt and Martin Nilsson

SIAM J. Sci. Comput. 25, pp. 2126-2144 (19 pages) | Cited 5 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
An efficient method for the solution of systems of linear equations with many right-hand sides is developed. The right-hand sides are assumed to depend smoothly on a parameter. The equations are solved by an iterative method, and a linear least squares approximation is used as initial guess. The work spent on the iterations is bounded independently of the number of right-hand sides. The method is applied to the solution of Maxwell's equations of electromagnetics in the frequency domain. The efficiency of the method is illustrated by computing the monostatic radar cross section around an aircraft model.

Long-Time-Step Integrators for Almost-Adiabatic Quantum Dynamics

Tobias Jahnke

SIAM J. Sci. Comput. 25, pp. 2145-2164 (20 pages) | Cited 4 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
The highly oscillatory solution of a singularly perturbed Schrödinger equation with time-dependent Hamiltonian is computed numerically. The new time-symmetric integrators presented here can be used efficiently with step sizes significantly larger than those required by traditional schemes. This is achieved by a transformation of the problem and an expansion technique for integrals over the oscillating components. The error behavior in the adiabatic case is thoroughly analyzed, and the performance of the methods is illustrated both in an almost-adiabatic setup and in an avoided energy level crossing, where nonadiabatic state transitions occur.

An HLLC-Type Approximate Riemann Solver for Ideal Magnetohydrodynamics

K. F. Gurski

SIAM J. Sci. Comput. 25, pp. 2165-2187 (23 pages) | Cited 18 times

Online Publication Date: July 25, 2006

Full Text: | Download PDF

Show Abstract
This paper presents a solver based on the HLLC (Harten--Lax--van Leer contact wave) approximate nonlinear Riemann solver for gas dynamics for the ideal magnetohydrodynamics (MHD) equations written in conservation form. It is shown how this solver also can be considered a modification of Linde's "adequate" solver. This approximation method is intended to resolve slow, Alfvén, and contact waves better than the original HLL (Harten--Lax--van Leer) solver. Compared to exact nonlinear solvers and Roe's solver, this new solver is computationally inexpensive. In addition, the method will exactly resolve isolated contacts and fast shocks. The method also preserves positive density and pressure with two caveats: first, the numerical signal velocities (the eigenvalues of the Roe average matrix) do not underestimate the physical signal velocities, and second, in a very few cases it may be required to change the wavespeeds of the Riemann fan for the underlying HLL method to guarantee positive pressures. These conditions are less restrictive on the definitions of the wavespeeds than the conditions needed to make the HLLC method positively conservative for gas dynamics. While the method is intended for a three-dimensional MHD problem, the simulation results concentrate on one-dimensional test cases.
Close

close