Showing papers in "SIAM Journal on Scientific Computing in 2004"

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

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.

Computing the Ground State Solution of Bose--Einstein Condensates by a Normalized Gradient Flow

Abstract: In this paper, we present a continuous normalized gradient flow (CNGF) and prove its energy diminishing property, which provides a mathematical justification of the imaginary time method used in the physics literature to compute the ground state solution of Bose--Einstein condensates (BEC). We also investigate the energy diminishing property for the discretization of the CNGF. Two numerical methods are proposed for such discretizations: one is the backward Euler centered finite difference (BEFD) method, the other is an explicit time-splitting sine-spectral (TSSP) method. Energy diminishing for BEFD and TSSP for the linear case and monotonicity for BEFD for both linear and nonlinear cases are proven. Comparison between the two methods and existing methods, e.g., Crank--Nicolson finite difference (CNFD) or forward Euler finite difference (FEFD), shows that BEFD and TSSP are much better in terms of preserving the energy diminishing property of the CNGF. Numerical results in one, two, and three dimensions with magnetic trap confinement potential, as well as a potential of a stirrer corresponding to a far-blue detuned Gaussian laser beam, are reported to demonstrate the effectiveness of BEFD and TSSP methods. Furthermore we observe that the CNGF and its BEFD discretization can also be applied directly to compute the first excited state solution in BEC when the initial data is chosen as an odd function.

An Extension of MATLAB to Continuous Functions and Operators

Abstract: An object-oriented MATLAB system is described for performing numerical linear algebra on continuous functions and operators rather than the usual discrete vectors and matrices. About eighty MATLAB functions from plot and sum to svd and cond have been overloaded so that one can work with our "chebfun" objects using almost exactly the usual MATLAB syntax. All functions live on [-1,1] and are represented by values at sufficiently many Chebyshev points for the polynomial interpolant to be accurate to close to machine precision. Each of our overloaded operations raises questions about the proper generalization of familiar notions to the continuous context and about appropriate methods of interpolation, differentiation, integration, zerofinding, or transforms. Applications in approximation theory and numerical analysis are explored, and possible extensions for more substantial problems of scientific computing are mentioned.

Numerical Simulation of the Smoluchowski Coagulation Equation

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.

An Automated Multilevel Substructuring Method for Eigenspace Computation in Linear Elastodynamics

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--Poincare 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.

Permuting Sparse Rectangular Matrices into Block-Diagonal Form

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.

Discrete Event Simulation of Hybrid Systems

Abstract: This paper describes the quantization-based integration methods and extends their use to the simulation of hybrid systems. Using the fact that these methods approximate ordinary differential equations (ODEs) and differential algebraic equations (DAEs) by discrete event systems, it is shown how hybrid systems can be approximated by pure discrete event simulation models (within the DEVS formalism framework). In this way, the treatment and detection of events representing discontinuities---which constitute an important problem for classic ODE solvers---is notably simplified. It can be also seen that the main advantages of quantization-based methods (error control, reduction of computational costs, possibilities of parallelization, sparsity exploitation, etc.) are still verified in the presence of discontinuities. Finally, some examples which illustrate the use and the advantages of the methodology in hybrid systems are discussed.

An HLLC-Type Approximate Riemann Solver for Ideal Magnetohydrodynamics

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, Alfven, 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 conditi...

Adaptive Smoothed Aggregation ($\alpha$SA)

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 Taxonomy of Consistently Stabilized Finite Element Methods for the Stokes Problem

Abstract: Stabilized mixed methods can circumvent the restrictive inf-sup condition without introducing penalty errors. For properly chosen stabilization parameters these methods are well-posed for all conforming velocity-pressure pairs. However, their variational forms have widely varying properties. First, stabilization offers a choice between weakly or strongly coercive bilinear forms that give rise to linear systems with identical solutions but very different matrix properties. Second, coercivity may be conditional upon a proper choice of a stabilizing parameter. Here we focus on how these two aspects of stabilized methods affect their accuracy and efficient iterative solution. We present results that indicate a preference of Krylov subspace solvers for strongly coercive formulations. Stability criteria obtained by finite element and algebraic analyses are compared with numerical experiments. While for two popular classes of stabilized methods, sufficient stability bounds correlate well with numerical stability, our experiments indicate the intriguing possibility that the pressure-stabilized Galerkin method is unconditionally stable.

Analysis of Preconditioners for Saddle-Point Problems

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.

Simulations of the Whirling Instability by the Immersed Boundary Method

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

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

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 objective is 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.

An Analysis of Delay-Dependent Stability for Ordinary and Partial Differential Equations with Fixed and Distributed Delays

Abstract: This paper is concerned with the study of the stability of ordinary and partial differential equations with both fixed and distributed delays, and with the study of the stability of discretizations of such differential equations. We start with a delay-dependent asymptotic stability analysis of scalar ordinary differential equations with real coefficients. We study the exact stability region of the continuous problem as a function of the parameters of the model. Next, it is proved that a time discretization based on the trapezium rule can preserve the asymptotic stability for the considered set of test problems. In the second part of the paper, we study delay partial differential equations. The stability region of the fully continuous problem is analyzed first. Then a semidiscretization in space is applied. It is shown that the spatial discretization leads to a reduction of the stability region when the standard second-order central difference operator is employed to approximate the diffusion operator. Finally we consider the delay-dependent stability of the fully discrete problem, where the partial differential equation is discretized both in space and in time. Some numerical examples and further discussions are given.

The Ultra-Weak Variational Formulation for Elastic Wave Problems

Abstract: The ultra-weak variational formulation has been used effectively to solve time-harmonic acoustic and electromagnetic wave propagation in inhomogeneous media. We develop the ultra-weak variational formulation for elastic wave propagation in two space dimensions. In order to improve the accuracy and stability of the method, we find it necessary to approximate the S- and P-wave components of the solution in a balanced way. Some preliminary analysis is provided and numerical evidence is presented for the efficiency of the scheme in comparison to piecewise linear finite elements.

An Implicit-Explicit Runge--Kutta--Chebyshev Scheme for Diffusion-Reaction Equations

Abstract: An implicit-explicit (IMEX) extension of the explicit Runge--Kutta--Chebyshev (RKC) scheme designed for parabolic PDEs is proposed for diffusion-reaction problems with severely stiff reaction terms. The IMEX scheme treats these reaction terms implicitly and diffusion terms explicitly. Within the setting of linear stability theory, the new IMEX scheme is unconditionally stable for reaction terms having a Jacobian matrix with a real spectrum. For diffusion terms the stability characteristics remain unchanged. A numerical comparison for a stiff, nonlinear radiation-diffusion problem between an RKC solver, an IMEX-RKC solver, and the popular implicit BDF solver VODPK using the Krylov solver GMRES illustrates the excellent performance of the new scheme.

Optimal Discrete Transmission Conditions for a Nonoverlapping Domain Decomposition Method for the Helmholtz Equation

Abstract: This paper is dedicated to recent developments of a two-Lagrange multipliers domain decomposition method for the Helmholtz equation [C. Farhat et al., Comput. Methods Appl. Mech. Engrg., 184 (2000), pp. 213--240; M. J. Gander, F. Magoules, and F. Nataf, SIAM J. Sci. Comput., 24 (2002), pp. 38--60] involving an additional augmented operator along the interface between the subdomains. Most methods for optimizing the augmented interface operator are based on the discretization of approximations of the continuous transparent operator [B. Despres, Proceedings of the Second International Conference on Mathematical and Numerical Aspects of Wave Propagation, R. Kleinman et al., eds., SIAM, Philadelphia, 1993, pp. 197--206; J.-D. Benamou and B. Despres, J. Comput. Phys., 136 (1997), pp. 68--82; P. Chevalier and F. Nataf, Domain Decomposition Methods 10, AMS, Providence, RI, 1998, pp. 400--407; M. J. Gander, Proceedings of the 12th International Conference on Domain Decomposition Methods, (Chiba, Japan), ddm.org, 2000, pp. 247--254; M. J. Gander, F. Magoules, and F. Nataf, SIAM J. Sci. Comput., 24 (2002), pp. 38--60]. At the discrete level, the optimal operator can be proved to be equal to the Schur complement of the outer domain. This Schur complement can be directly approximated using purely algebraic techniques like sparse approximate inverse methods or incomplete factorization. The main advantage of such an algebraic approach is that it is much easier to implement in existing code without any information on the geometry of the interface and the finite element formulation used. Convergence results and parallel efficiency of several algebraic optimization techniques of an interface operator for acoustic analysis applications will be presented.

Direct Evaluation of Hypersingular Galerkin Surface Integrals

Abstract: A direct algorithm for evaluating hypersingular integrals arising in a three-dimensional Galerkin boundary integral analysis is presented. The singular integrals are defined as limits to the boundary, and by integrating two of the four dimensions analytically, the coincident integral is shown to be divergent. However, the divergent terms can be explicitly calculated and shown to cancel with corresponding singularities in the adjacent edge integrals. A single analytic integration is employed for the edge and vertex singular integrals. This is sufficient to display the divergent term in the edge-adjacent integral and to show that the vertex integral is finite. By explicitly identifying the divergent quantities, we can compute the hypersingular integral without recourse to Stokes's theorem or the Hadamard finite part. The algorithms are developed in the context of a linear element approximation for the Laplace equation but are expected to be generally applicable. As an example, the algorithms are applied to solve a thermal problem in an exponentially graded material.

A Jacobi--Davidson Method for Solving Complex Symmetric Eigenvalue Problems

Abstract: We discuss variants of the Jacobi--Davidson method for solving the generalized complex symmetric eigenvalue problem. The Jacobi--Davidson algorithm can be considered as an accelerated inexact Rayleigh quotient iteration. We show that it is appropriate to replace the Euclidean inner product in ${\mathbb C}^n$ with an indefinite inner product. The Rayleigh quotient based on this indefinite inner product leads to an asymptotically cubically convergent Rayleigh quotient iteration. Advantages of the method are illustrated by numerical examples. We deal with problems from electromagnetics that require the computation of interior eigenvalues. The main drawback that we experience in these particular examples is the lack of efficient preconditioners.

A Jump Condition Capturing Finite Difference Scheme for Elliptic Interface Problems

Abstract: We propose a simple finite difference scheme for the elliptic interface problem with a discontinuous diffusion coefficient using a body-fitted curvilinear coordinate system The resulting matrix is symmetric and positive definite Standard techniques of acceleration such as PCG and multigrid can be used to invert the matrix The main advantage of the scheme is its simplicity: the entries of the matrix are simply the centered difference second order approximation of the metric tensor $g^{\alpha\beta}$ In addition, the interface jump conditions are naturally built into the finite difference discretization No interpolation/extrapolation process is involved in the derivation of the scheme Both the solution and the flux are observed to have second order accuracy

Partially Stirred Reactor Model: Analytical Solutions and Numerical Convergence Study of a PDF/Monte Carlo Method

Abstract: We investigate the partially stirred reactor (PaSR), which is based on a simplified joint composition probability density function (PDF) transport equation Analytical solutions for the first four moments of the mass density function (MDF) obtained from the PaSR model are presented The Monte Carlo particle method with first order time splitting algorithm is implemented to obtain the first four moments of the MDF numerically The dynamics of the stochastic particle system is determined by inflow-outflow, chemical reaction, and mixing events Three different inflow-outflow algorithms are investigated: an algorithm based on the inflow-outflow event modeled as a Poisson process, an inflow-outflow algorithm mentioned in the literature, and a novel algorithm derived on the basis of analytical solutions It is demonstrated that the inflow-outflow algorithm used in the literature can be explained by considering a deterministic waiting time parameter of a corresponding stochastic process, and also forms a specific case of the new algorithm The number of particles in the ensemble, N, the nondimensional time step, $\Delta t^{*}$ (ratio of the global time step to the characteristic time of an event), and the number of independent simulation trials, L, are the three sources of the numerical error The split analytical solutions and the numerical experiments suggest that the systematic error converges as $\Delta t^{*}$ and N-1 The statistical error scales as L-1/2 and N-1/2 The significance of the numerical parameters and the inflow-outflow algorithms is also studied by applying the PaSR model to a practical case of premixed kerosene and air combustion

New Geometric Immersed Interface Multigrid Solvers

Abstract: In [L. Adams and Z. Li, SIAM J. Sci. Comput., 24 (2002), pp. 463--479], a multigrid method was designed specifically for interface problems that have been discretized using the methods described therein, and in [Z. Li and K. Ito, SIAM J. Sci. Comput., 23 (2001), pp. 339--361] for elliptic interface problems using the maximum principle preserving schemes. In this paper, we improve on this method by giving a new interpolator for grid points near the immersed interface and a new restrictor that guarantees the coarse-grid matrices are M-matrices. We compare this new restrictor to injection and the transpose of interpolation. We show that the number of V-cycles is constant as the mesh size decreases and increases only slightly as the ratio of the discontinuous problem coefficient grows at the interface only when this new restrictor is used.

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

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.

On Monotonicity of Difference Schemes for Computational Physics

Abstract: Criteria are developed for monotonicity of linear as well as nonlinear difference schemes associated with the numerical analysis of systems of partial differential equations, integrodifferential equations, etc. Difference schemes are converted into variational forms that satisfy the boundary maximum principle and also allow the investigation of monotonicity for nonlinear operators using linear patterns. Sufficient conditions are provided to review the monotonicity of single and coupled difference schemes. Necessary as well as necessary and sufficient conditions for monotonicity of explicit schemes are also developed. The notion of submonotone difference schemes is considered and the associated criteria are developed. We discuss the interrelationship between monotonicity, submonotonicity, and stability. Some known schemes serve as examples demonstrating the implementation of the developed approaches. Among these examples, we describe the possibility that stable schemes such as total variation diminishing (TVD) as well as monotonicity preserving can produce spurious oscillations.

Long-time-step integrators for almost-adiabatic quantum dynamics ∗

Abstract: The highly oscillatory solution of a singularly perturbed Schrodinger 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.

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

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

Linear versus Nonlinear Dimensionality Reduction of High-Dimensional Dynamical Systems

Abstract: Two techniques for dimensionality reduction of high-dimensional dynamical systems are presented. The first is based on Karhunen--Loeve (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 r...

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

Abstract: A minimum residual interpolation method for linear equations with multiple right hand sides

An Adaptive Mixed Scheme for Energy-Transport Simulations of Field-Effect Transistors

Abstract: Energy-transport models are used in semiconductor simulations to account for thermal effects. The model consists of the continuity equations for the number and energy of the electrons, coupled to the Poisson equation for the electrostatic potential. The movement of the holes is modeled by drift-diffusion equations, and Shockley--Read--Hall recombination-generation processes are taken into account. The stationary equations are discretized using a mixed-hybrid finite-element method introduced by Marini and Pietra. The two-dimensional mesh is adaptively refined using an error estimator motivated by results of Hoppe and Wohlmuth. The numerical scheme is applied to the simulation of a two-dimensional double-gate MESFET and a deep submicron MOSFET device.

Preconditioning Strategies for Hermitian Indefinite Toeplitz Linear Systems

Abstract: In this paper we propose and analyze preconditioning strategies for Hermitian indefinite linear systems by using indefinite preconditioners: under very elementary assumptions, we show that the eigenvalues are real. Moreover, in the case of multilevel Toeplitz structures, we prove distributional and localization results. These techniques used in connection with the CG, GMRES, BICGstab, and QMR algorithms allow us to solve in an optimal way the corresponding linear systems. A wide numerical experimentation confirms the efficiency of the proposed procedures.