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 Next Issue

2006

Volume 28, Issue 5, pp. 1597-1999


A Numerical Study of the Probe Method

Klaus Erhard and Roland Potthast

SIAM J. Sci. Comput. 28, pp. 1597-1612 (16 pages) | Cited 5 times

Online Publication Date: September 29, 2006

Full Text: | Download PDF

Show Abstract
The goal of this work is the numerical realization of the probe method suggested by Ikehata for the detection of an obstacle D in inverse scattering. The main idea of the method is to use probes in the form of point source $\Phi(\cdot,z)$ with source point z to define an indicator function $\hat{I}(z)$ which can be reconstructed from Cauchy data or far field data. The indicator function $\hat{I}(z)$ can be shown to blow off when the source point z tends to the boundary $\partial D$, and this behavior can be used to find D. To study the feasibility of the probe method we will use two equivalent formulations of the indicator function. We will carry out the numerical realization of the functional and show reconstructions of a sound‐soft obstacle.

Benefits of IEEE‐754 Features in Modern Symmetric Tridiagonal Eigensolvers

Osni A. Marques, E. Jason Riedy, and Christof Vömel

SIAM J. Sci. Comput. 28, pp. 1613-1633 (21 pages) | Cited 1 time

Online Publication Date: September 29, 2006

Full Text: | Download PDF

Show Abstract
Bisection is one of the most common methods used to compute the eigenvalues of symmetric tridiagonal matrices. Bisection relies on the Sturm count: For a given shift σ, the number of negative pivots in the factorization $T - \sigma I = LDL^T$ equals the number of eigenvalues of T that are smaller than σ. In IEEE‐754 arithmetic, the value ∞ permits the computation to continue past a zero pivot, producing a correct Sturm count when T is unreduced. Demmel and Li showed [IEEE Trans. Comput., 43 (1994), pp. 983–992] that using ∞ rather than testing for zero pivots within the loop could significantly improve performance on certain architectures. When eigenvalues are to be computed to high relative accuracy, it is often preferable to work with $LDL^T$ factorizations instead of the original tridiagonal T. One important example is the MRRR algorithm. When bisection is applied to the factored matrix, the Sturm count is computed from $LDL^T$ which makes differential stationary and progressive qds algorithms the methods of choice. While it seems trivial to replace T by $LDL^T$, in reality these algorithms are more complicated: In IEEE‐754 arithmetic, a zero pivot produces an overflow followed by an invalid exception (NaN, or “Not a Number”) that renders the Sturm count incorrect. We present alternative, safe formulations that are guaranteed to produce the correct result. Benchmarking these algorithms on a variety of platforms shows that the original formulation without tests is always faster provided that no exception occurs. The transforms see speed‐ups of up to $2.6\times$ over the careful formulations. Tests on industrial matrices show that encountering exceptions in practice is rare. This leads to the following design: First, compute the Sturm count by the fast but unsafe algorithm. Then, if an exception occurs, recompute the count by a safe, slower alternative. The new Sturm count algorithms improve the speed of bisection by up to $2\times$ on our test matrices. Furthermore, unlike the traditional tiny‐pivot substitution, proper use of IEEE‐754 features provides a careful formulation that imposes no input range restrictions.

A Cache‐Aware Algorithm for PDEs on Hierarchical Data Structures Based on Space‐Filling Curves

Frank Günther, Miriam Mehl, Markus Pögl, and Christoph Zenger

SIAM J. Sci. Comput. 28, pp. 1634-1650 (17 pages) | Cited 6 times

Online Publication Date: October 06, 2006

Full Text: | Download PDF

Show Abstract
Competitive numerical algorithms for solving partial differential equations have to work with the most efficient numerical methods like multigrid and adaptive grid refinement and thus with hierarchical data structures. Unfortunately, in most implementations, hierarchical data—typically stored in trees—cause a nonnegligible overhead in data access. To overcome this quandary—numerical efficiency versus efficient implementation—our algorithm uses space‐filling curves to build up data structures which are processed linearly. In fact, the only kind of data structure used in our implementation is stacks. Thus, data access becomes very fast—even faster than the common access to nonhierarchical data stored in matrices—and, in particular, cache misses are reduced considerably. Furthermore, the implementation of multigrid cycles and/or higher order discretizations as well as the parallelization of the whole algorithm become very easy and straightforward on these data structures.

Recycling Krylov Subspaces for Sequences of Linear Systems

Michael L. Parks, Eric de Sturler, Greg Mackey, Duane D. Johnson, and Spandan Maiti

SIAM J. Sci. Comput. 28, pp. 1651-1674 (24 pages) | Cited 29 times

Online Publication Date: October 06, 2006

Full Text: | Download PDF

Show Abstract
Many problems in science and engineering require the solution of a long sequence of slowly changing linear systems. We propose and analyze two methods that significantly reduce the total number of matrix‐vector products required to solve all systems. We consider the general case where both the matrix and right‐hand side change, and we make no assumptions regarding the change in the right‐hand sides. Furthermore, we consider general nonsingular matrices, and we do not assume that all matrices are pairwise close or that the sequence of matrices converges to a particular matrix. Our methods work well under these general assumptions, and hence form a significant advancement with respect to related work in this area. We can reduce the cost of solving subsequent systems in the sequence by recycling selected subspaces generated for previous systems. We consider two approaches that allow for the continuous improvement of the recycled subspace at low cost. We consider both Hermitian and non‐Hermitian problems, and we analyze our algorithms both theoretically and numerically to illustrate the effects of subspace recycling. We also demonstrate the effectiveness of our algorithms for a range of applications from computational mechanics, materials science, and computational physics.

On Mass‐Conserving Least‐Squares Methods

J. J. Heys, E. Lee, T. A. Manteuffel, and S. F. Mccormick

SIAM J. Sci. Comput. 28, pp. 1675-1693 (19 pages) | Cited 6 times

Online Publication Date: October 06, 2006

Full Text: | Download PDF

Show Abstract
Least‐squares variational methods have several practical and theoretical advantages for solving elliptic partial differential equations, including symmetric positive definite discrete operators and a sharp error measure. One of the potential drawbacks, especially in three dimensions, is that mass conservation is achieved only in a least‐squares sense, and underresolved solutions are especially vulnerable to poor conservation. For the stationary Navier–Stokes equations, which are typically rewritten as a larger system of first‐order equations, the loss of mass in the approximate solution is strongly dependent upon the boundary conditions used. A new first‐order system reformulation of the Navier–Stokes equations is presented that admits a wider range of mass‐conserving boundary conditions. This new formulation is shown to provide both excellent mass conservation and excellent algebraic multigrid performance for three different problems, a square channel, backward‐facing step, and branching tubes with two generations of bifurcations.

A Variational Approach to Modeling and Simulation of Grain Growth

David Kinderlehrer, Irene Livshits, and Shlomo Ta’asan

SIAM J. Sci. Comput. 28, pp. 1694-1715 (22 pages) | Cited 5 times

Online Publication Date: October 12, 2006

Full Text: | Download PDF

Show Abstract
Most technologically useful materials arise as polycrystalline microstructures, composed of a myriad of small crystallites, called grains, separated by their interfaces, called grain boundaries. The orientations and arrangements of the grains and their network of boundaries are implicated in many properties across wide scales, for example, functional properties, like conductivity in microprocessors, and lifetime properties, like fracture toughness in structures. Simulation is becoming an important tool for understanding both materials properties and their processing requirements. Here we offer a consistent variational approach to the mesoscale simulation of these systems subject to the Mullins equation of curvature‐driven growth in a two‐dimensional setting. The main objective is to provide a calibration for future two‐dimensional and three‐dimensional efforts. We discuss several novel features of our approach, which we anticipate will render it a flexible, scalable, and robust tool to aid in microstructural prediction. Simulation results offer compelling evidence of the predictability and robustness of statistical properties of large systems, such as grain size distribution and texture, that are of immediate interest in materials science.

A Pseudospectral Fictitious Point Method for High Order Initial‐Boundary Value Problems

Bengt Fornberg

SIAM J. Sci. Comput. 28, pp. 1716-1729 (14 pages) | Cited 3 times

Online Publication Date: October 16, 2006

Full Text: | Download PDF

Show Abstract
When pseudospectral approximations are used for space derivatives, one often encounters spurious eigenvalues. These can lead to severe time stepping difficulties for PDEs. This is especially the case for equations with high order derivatives in space, requiring multiple conditions at one or both boundaries. We note here that a very simple‐to‐implement fictitious point approach circumvents most of these difficulties. The new approach is tested on the Kuramoto–Sivashinsky equation and on a dispersive linear PDE featuring a time‐space corner singularity.

A High‐Order Accurate Parallel Solver for Maxwell’s Equations on Overlapping Grids

William D. Henshaw

SIAM J. Sci. Comput. 28, pp. 1730-1765 (36 pages) | Cited 3 times

Online Publication Date: October 20, 2006

Full Text: | Download PDF

Show Abstract
A scheme for the solution of the time‐dependent Maxwell’s equations on composite overlapping grids is described. The method uses high‐order accurate approximations in space and time for Maxwell’s equations written as a second‐order vector wave equation. High‐order accurate symmetric difference approximations to the generalized Laplace operator are constructed for curvilinear component grids. The modified equation approach is used to develop high‐order accurate approximations that use only three time levels and have the same time‐stepping restriction as the second‐order scheme. Discrete boundary conditions for perfect electrical conductors and for material interfaces are developed and analyzed. The implementation is optimized for component grids that are Cartesian, resulting in a fast and efficient method. The solver runs on parallel machines with each component grid distributed across one or more processors. Numerical results in two and three dimensions are presented for the fourth‐order accurate version of the method. These results demonstrate the accuracy and efficiency of the approach.

An Unstaggered, High‐Resolution Constrained Transport Method for Magnetohydrodynamic Flows

James A. Rossmanith

SIAM J. Sci. Comput. 28, pp. 1766-1797 (32 pages) | Cited 10 times

Online Publication Date: October 20, 2006

Full Text: | Download PDF

Show Abstract
The ideal magnetohydrodynamic (MHD) equations are important in modeling phenomena in a wide range of applications, including space weather, solar physics, laboratory plasmas, and astrophysical fluid flows. Numerical methods for the MHD equations must confront the challenge of producing approximate solutions that remain accurate near shock waves and that satisfy a divergence‐free constraint on the magnetic field. Failure to accomplish this often leads to unphysical solutions. In this paper, a high‐resolution wave propagation method is developed that utilizes a novel constrained transport technique to keep the magnetic field divergence‐free. This approach is based on directly solving the magnetic potential equation in conjunction with a new limiting strategy to obtain a nonoscillatory magnetic field. It is demonstrated in this work that an unstaggered definition of the divergence is the correct one to use in the case of wave propagation methods. Therefore, we solve the magnetic potential equation on the same grid as the MHD equations; hence the usual grid staggering that is found in constrained transport methods is eliminated. We demonstrate through truncation error analysis and direct numerical simulation that the resulting method is second order accurate in space and time for smooth solutions and nonoscillatory near shocks and other discontinuities. The resulting numerical method has been implemented as an extension to the clawpack software package and can be freely downloaded from the Web.

A Rational Spectral Collocation Method with Adaptively Transformed Chebyshev Grid Points

T. W. Tee and Lloyd N. Trefethen

SIAM J. Sci. Comput. 28, pp. 1798-1811 (14 pages) | Cited 25 times

Online Publication Date: October 24, 2006

Full Text: | Download PDF

Show Abstract
A spectral collocation method based on rational interpolants and adaptive grid points is presented. The rational interpolants approximate analytic functions with exponential accuracy by using prescribed barycentric weights and transformed Chebyshev points. The locations of the grid points are adapted to singularities of the underlying solution, and the locations of these singularities are approximated by the locations of poles of Chebyshev–Padé approximants. Numerical experiments on two time‐dependent problems, one with finite time blow‐up and one with a moving front, indicate that the method far outperforms the standard Chebyshev spectral collocation method for problems whose solutions have singularities in the complex plane close to $[-1,1]$.

Efficient Computation of Isometry‐Invariant Distances Between Surfaces

Alexander M. Bronstein, Michael M. Bronstein, and Ron Kimmel

SIAM J. Sci. Comput. 28, pp. 1812-1836 (25 pages) | Cited 27 times

Online Publication Date: October 27, 2006

Full Text: | Download PDF

Show Abstract
We present an efficient computational framework for isometry‐invariant comparison of smooth surfaces. We formulate the Gromov–Hausdorff distance as a multidimensional scaling–like continuous optimization problem. In order to construct an efficient optimization scheme, we develop a numerical tool for interpolating geodesic distances on a sampled surface from precomputed geodesic distances between the samples. For isometry‐invariant comparison of surfaces in the case of partially missing data, we present the partial embedding distance, which is computed using a similar scheme. The main idea is finding a minimum‐distortion mapping from one surface to another, while considering only relevant geodesic distances. We discuss numerical implementation issues and present experimental results that demonstrate its accuracy and efficiency.

Extensions of Priest’s Double‐Precision Summation

Yves Nievergelt

SIAM J. Sci. Comput. 28, pp. 1837-1850 (14 pages)

Online Publication Date: November 14, 2006

Full Text: | Download PDF

Show Abstract
From single‐precision summands sorted in nonincreasing exponents, extended‐ precision faithful additions with any radix yield an extended‐precision sum, which, rounded back to single‐precision, suffers from a roundoff smaller than one unit in the penultimate significant digit. Such provable accuracy holds only if the radix‐based logarithm of the number of summands does not exceed the difference between the number of digits in extended and single precisions.

Wavelet Filtering for Exact Controllability of the Wave Equation

M. Negreanu, A.‐M. Matache, and C. Schwab

SIAM J. Sci. Comput. 28, pp. 1851-1885 (35 pages) | Cited 3 times

Online Publication Date: November 14, 2006

Full Text: | Download PDF

Show Abstract
In this paper we propose and analyze a wavelet‐based numerical method for exact controllability of the one‐dimensional wave equation, concentrating on the particular case of Dirichlet controls. We use the Hilbert uniqueness method (HUM) and prove this result in the context of the finite element space semidiscretization with spline wavelet basis functions. The algorithm for the numerical construction of the controls is based on a conjugate gradient method proposed by Glowinski [J. Comput. Phys., 103 (1992), pp. 189–221] and a two‐grid filtering technique in the wavelet basis. We show that wavelet filtering yields numerically exact controllability in the optimal time of the continuous problem and the convergence of the numerical controls using a multilevel discretization filtering technique. Classical arguments then allow proving the uniform boundedness of the numerical controls and passing to the limit as the mesh size tends to zero. We prove that the condition number of the discretized HUM operator $\Lambda$ with wavelet filtering is uniformly bounded with respect to the mesh width and prove convergence of the numerical controls obtained with the two‐grid filtering technique. Numerical experiments show that wavelet filtering yields numerical controls with the same control time $T$ as in the continuous case.

A Parallel Implementation of Dual‐Primal FETI Methods for Three‐Dimensional Linear Elasticity Using a Transformation of Basis

Axel Klawonn and Oliver Rheinbach

SIAM J. Sci. Comput. 28, pp. 1886-1906 (21 pages) | Cited 15 times

Online Publication Date: November 16, 2006

Full Text: | Download PDF

Show Abstract
Dual‐primal FETI methods for linear elasticity problems in three dimensions are considered. These are nonoverlapping domain decomposition methods where some primal continuity constraints across subdomain boundaries are required to hold throughout the iterations, whereas most of the constraints are enforced by Lagrange multipliers. An algorithmic framework for dual‐primal FETI methods is described together with a transformation of basis to implement the primal constraints. Numerical results obtained from a parallel implementation of these algorithms applied to a model benchmark problem with structured meshes and to problems with more complicated geometries from industrial and biological applications using unstructured meshes are provided. These results show that the presented dual‐primal FETI algorithms are numerical and parallel scalable.

Parallel Guaranteed Quality Delaunay Uniform Mesh Refinement

Andrey N. Chernikov and Nikos P. Chrisochoides

SIAM J. Sci. Comput. 28, pp. 1907-1926 (20 pages) | Cited 1 time

Online Publication Date: November 16, 2006

Full Text: | Download PDF

Show Abstract
We present a theoretical framework for developing parallel guaranteed quality Delaunay mesh generation software that allows us to use commercial off‐the‐shelf sequential Delaunay meshers for two‐dimensional geometries. In this paper, we describe our approach for constructing uniform meshes, that is, the meshes in which all elements have approximately the same size. Our uniform distributed‐ and shared‐memory implementations are based on a simple (block) coarse‐grained mesh decomposition. Our method requires only local communication, which is bulk and structured as opposed to fine and unpredictable communication of the other existing practical parallel guaranteed quality mesh generation and refinement techniques. Our experimental data show that on a cluster of more than 100 workstations we can generate about 0.9 billion elements in less than 5 minutes in the absence of work‐load imbalances. Preliminary results for this paper were presented in [A. N. Chernikov and N. P. Chrisochoides, “Practical and efficient point insertion scheduling method for parallel guaranteed quality Delaunay refinement,” in Proceedings of the 18th Annual International Conference on Supercomputing, ACM Press, New York, 2004, pp. 48–57]. Our work in progress includes extending the presented approach, which can efficiently generate only uniform meshes, to nonuniform graded meshes.

Behavior of Finite Volume Schemes for Hyperbolic Conservation Laws on Adaptive Redistributed Spatial Grids

Ch. Arvanitis and A. I. Delis

SIAM J. Sci. Comput. 28, pp. 1927-1956 (30 pages)

Online Publication Date: November 16, 2006

Full Text: | Download PDF

Show Abstract
In this work we consider finite volume schemes combined with dynamic spatial mesh redistribution. We study whether appropriate mesh redistribution is a satisfactory mechanism for increasing the resolution of numerical solutions for problems of scalar and systems of conservation laws (CL) in one space dimension, while being at the same time a stabilization mechanism for selecting the appropriate entropy solution. In order to increase the resolution around shock areas and keep the computational cost low, our redistribution policy is to reconstruct spatially the numerical solution on a new mesh, where the solution’s curvature is almost uniformly distributed, while the node’s cardinality is kept constant. We examine the stabilization properties of that redistribution process by adding it as a substep on the time evolution step of some classical schemes with known (unstable) characteristics. Testing the resulting method for several such schemes and on a large number of CL problems that have solutions with special characteristics (shocks, rarefaction areas, steady states) and comparing the results with those produced by schemes with extra stabilization mechanisms (like slope/flux limiters, entropy corrections), we conclude that indeed the proposed redistribution adds such stabilization properties while at the same time increasing the resolution.

Segmentation with Depth: A Level Set Approach

Wei Zhu, Tony Chan, and Selim Esedoḡlu

SIAM J. Sci. Comput. 28, pp. 1957-1973 (17 pages) | Cited 6 times

Online Publication Date: November 24, 2006

Full Text: | Download PDF

Show Abstract
Segmentation with depth is the challenging problem of obtaining three dimensional information from a single two dimensional image. Unlike the standard segmentation problem, the goal of segmentation with depth is to determine not only the boundaries of objects that appear in the image, but also their relative distances to the observer by making use of occlusion relations. Nitzberg, Mumford, and Shiota [The 2.1D sketch, in Proceedings of the Third International Conference on Computer Vision, Osaka, Japan, 1990; Filtering, Segmentation, and Depth, Lecture Notes in Comput. Sci. 662, Springer‐Verlag, Berlin, 1993] proposed a variational formulation of this problem; according to their model, called the 2.1D sketch model, the regions that the objects occupy and their relative depth are to be extracted from the two dimensional image by minimizing a curvature based functional. Numerically, this is a highly nontrivial problem as the functional involves curvatures of the unknown contours. In this paper, we develop a level set based procedure for minimizing the Nitzberg–Mumford–Shiota energy. Unlike the minimization technique of Nitzberg–Mumford–Shiota [Filtering, Segmentation, and Depth, Lecture Notes in Comput. Sci. 662, Springer‐Verlag, Berlin, 1993], our technique represents contours implicitly and thus allows for connections between T‐junctions to take place automatically.

On the Acoustic Single Layer Potential: Stabilization and Fourier Analysis

A. Buffa and S. Sauter

SIAM J. Sci. Comput. 28, pp. 1974-1999 (26 pages) | Cited 4 times

Online Publication Date: November 24, 2006

Full Text: | Download PDF

Show Abstract
In this paper, we propose a general approach for stabilizing the single layer potential for the Helmholtz boundary integral equation and prove its stability. We consider Galerkin boundary element discretizations and analyze their convergence. Furthermore, we derive quantitative error bounds for the Galerkin discretization which are explicit with respect to the mesh width and the wave number for the special case that the surface is the unit sphere in $% \mathbb{R}^3$. We perform then a qualitative analysis which allows us to choose the stabilization such that the (negative) influence of the wave number in the stability and convergence estimates attains its minumum.
Close

Your Recent History Recent History Help

Recently Viewed

  1. A New Error Estimate of the Fast Gauss Transform
  2. A Preconditioner for the Steady-State Navier--Stokes Equations
  3. Iterative Solution of Cyclically Reduced Systems Arising from Discretization of the Three-Dimensiona...
  4. Nonoscillatory Central Schemes for Multidimensional Hyperbolic Conservation Laws
  5. A Fast Multilevel Algorithm for the Solution of Nonlinear Systems of Conductive-Radiative Heat Trans...
© 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