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

2008

Volume 31, Issue 2, pp. 799-1602


A Discrete Boltzmann Equation Based on a Cub-Octahedron in $\mathbb{R}^3$

Laek S. Andallah and Hans Babovsky

SIAM J. Sci. Comput. 31, pp. 799-825 (27 pages)

Online Publication Date: November 14, 2008

Full Text: | Download PDF

Show Abstract
In this article we introduce a discrete Boltzmann equation in $\mathbb{R}^3$ velocity space based on a cub-octahedron. We present a layer-wise construction of a regular collision model in a bounded grid of the sphere packing problem. The paper contains the analytical results about the basic kinetic features of the regular collision model. We analyze the computational complexity, report the efficiency of the model, and perform numerical experiments based on the discrete model Boltzmann equation for some standard test problem.

A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions

F. Fang and C. W. Oosterlee

SIAM J. Sci. Comput. 31, pp. 826-848 (23 pages) | Cited 2 times

Online Publication Date: November 14, 2008

Full Text: | Download PDF

Show Abstract
Here we develop an option pricing method for European options based on the Fourier-cosine series and call it the COS method. The key insight is in the close relation of the characteristic function with the series coefficients of the Fourier-cosine expansion of the density function. In most cases, the convergence rate of the COS method is exponential and the computational complexity is linear. Its range of application covers underlying asset processes for which the characteristic function is known and various types of option contracts. We will present the method and its applications in two separate parts. The first one is this paper, where we deal with European options in particular. In a follow-up paper we will present its application to options with early-exercise features.

Automated Code Generation for Discontinuous Galerkin Methods

Kristian B. Ølgaard, Anders Logg, and Garth N. Wells

SIAM J. Sci. Comput. 31, pp. 849-864 (16 pages) | Cited 3 times

Online Publication Date: November 21, 2008

Full Text: | Download PDF

Show Abstract
A compiler approach for generating low-level computer code from high-level input for discontinuous Galerkin finite element forms is presented. The input language mirrors conventional mathematical notation, and the compiler generates efficient code in a standard programming language. This facilitates the rapid generation of efficient code for general equations in varying spatial dimensions. Key concepts underlying the compiler approach and the automated generation of computer code are elaborated. The approach is demonstrated for a range of common problems, including the Poisson, biharmonic, advection-diffusion, and Stokes equations.

A Newton–Krylov Solver for Remapping-Based Volume-of-Fluid Methods

Petar Liovic and Djamel Lakehal

SIAM J. Sci. Comput. 31, pp. 865-889 (25 pages)

Online Publication Date: November 21, 2008

Full Text: | Download PDF

Show Abstract
A new conservative remapping method is presented here for generating multi-dimensional fluxes used for the time integration of hyperbolic equations through single-stage unsplit advection. A solution to the Lagrangian mesh location is generated that converges the volumes of all remapped flux definitions to the volumes of equivalent one-dimensional flux definitions at their corresponding cell faces. The method is most feasible within localized subdomains of an Eulerian solution domain. The main utility of the method is the generation of fluxes in the vicinity of discontinuities and sharp gradients that does not introduce physically spuriously extrema while locally conserving mass. In generating solutions to Lagrangian mesh location by solving a nonlinear equation system, the new remap requires the anchoring of the Lagrangian mesh—removal of degrees of freedom to make the problem definite. The method is based on characteristic tracing along trajectories assumed to be linear. This paper investigates the convergence and sensitivity of the Lagrangian mesh solver to mesh and timestep size, in assessing its merit for application in transient CFD codes. For single-stage unsplit advection volume-of-fluid methods, the new remap suppresses the introduction of undershoots, overshoots, and wisps. In combination with a conservative remapping phase, the new method for fixing the location of the Lagrangian mesh eliminates the need for heuristic redistribution or clipping procedures for spurious extrema removal.

Probing the Pareto Frontier for Basis Pursuit Solutions

Ewout van den Berg and Michael P. Friedlander

SIAM J. Sci. Comput. 31, pp. 890-912 (23 pages) | Cited 26 times

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
The basis pursuit problem seeks a minimum one-norm solution of an underdetermined least-squares problem. Basis pursuit denoise (BPDN) fits the least-squares problem only approximately, and a single parameter determines a curve that traces the optimal trade-off between the least-squares fit and the one-norm of the solution. We prove that this curve is convex and continuously differentiable over all points of interest, and show that it gives an explicit relationship to two other optimization problems closely related to BPDN. We describe a root-finding algorithm for finding arbitrary points on this curve; the algorithm is suitable for problems that are large scale and for those that are in the complex domain. At each iteration, a spectral gradient-projection method approximately minimizes a least-squares problem with an explicit one-norm constraint. Only matrix-vector operations are required. The primal-dual solution of this problem gives function and derivative information needed for the root-finding method. Numerical experiments on a comprehensive set of test problems demonstrate that the method scales well to large problems.

Discontinuous-Continuous Galerkin Methods for a Structured Model of a Biological System

Mi-Young Kim

SIAM J. Sci. Comput. 31, pp. 913-938 (26 pages)

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
A hybrid discontinuous-continuous Galerkin finite element method is proposed to approximate the solution to a class of composite hyperbolic-parabolic PDEs which arises in the study of biological systems. We employ a discontinuous Galerkin finite element technique for the numerical approximation of the linear hyperbolic transport part of the PDE posed in an age-time domain, while we approximate the second-order elliptic operator using a standard conforming finite element method. The strong stability is established, and a priori $L^2(L^2)$-error estimates are obtained. The finite elements in the age-time domain are constructed in such a manner that the discontinuous Galerkin method is applied over a triangulation in an explicit fashion from triangle to triangle. The methods are effective in the sense that age and time steppings are easily adaptive and the computations are readily parallelizable. No restriction is imposed on the time step in connection with the mesh size in the space. The approximate solution is computed slap by slap marching in time. Some numerical examples are presented.

Inertia-Revealing Preconditioning For Large-Scale Nonconvex Constrained Optimization

Olaf Schenk, Andreas Wächter, and Martin Weiser

SIAM J. Sci. Comput. 31, pp. 939-960 (22 pages) | Cited 1 time

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
Fast nonlinear programming methods following the all-at-once approach usually employ Newton's method for solving linearized Karush–Kuhn–Tucker (KKT) systems. In nonconvex problems, the Newton direction is guaranteed to be a descent direction only if the Hessian of the Lagrange function is positive definite on the nullspace of the active constraints; otherwise some modifications to Newton's method are necessary. This condition can be verified using the signs of the KKT eigenvalues (inertia), which are usually available from direct solvers for the arising linear saddle point problems. Iterative solvers are mandatory for very large-scale problems, but in general they do not provide the inertia. Here we present a preconditioner based on a multilevel incomplete $LBL^T$ factorization, from which an approximation of the inertia can be obtained. The suitability of the heuristics for application in optimization methods is verified on an interior point method applied to the CUTE and COPS test problems, on large-scale three-dimensional (3D) PDE-constrained optimal control problems, and on 3D PDE-constrained optimization in biomedical cancer hyperthermia treatment planning. The efficiency of the preconditioner is demonstrated on convex and nonconvex problems with $150^3$ state variables and $150^2$ control variables, both subject to bound constraints.

Stability Analysis of Time Stepping for Prolonged Plasma Fluid Simulations

L. Stals, R. Numata, and R. Ball

SIAM J. Sci. Comput. 31, pp. 961-986 (26 pages)

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
The Hasegawa–Wakatani system of equations may be used to predict and study the behavior of plasma flow. Stability analysis of the flow requires results over prolonged time series, which places a great strain on computational resources. Results can only be achieved for a wide choice of parameters by using numerical methods that allow long time steps and do not pollute the results with numerical instabilities. The report presents an analysis of several linear multistep methods and concludes that much of the understanding of the stability of linear systems also applies to the study of nonlinear problems such as the Hasegawa–Wakatani system of equations. In particular, methods such as the backward differentiation formulas should be used with the stiff systems generated by the discrete formulation of the Hasegawa–Wakatani system of equations.

Multibody Interactions in Coarse-Graining Schemes for Extended Systems

Sasanka Are, Markos A. Katsoulakis, Petr Plecháč, and Luc Rey-Bellet

SIAM J. Sci. Comput. 31, pp. 987-1015 (29 pages)

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
In this paper we address the role of multibody interactions for the coarse-grained approximation of stochastic lattice systems. Such interaction potentials are often not included in coarse-graining schemes, as they can be computationally expensive. The multibody interactions are obtained from the error expansion of the reference measure which is, in many cases, chosen as a Gibbs measure corresponding to a local mean-field approximation. We identify the parameter $\epsilon$ that characterizes the level of approximation and its relation to the underlying interaction potential. The error analysis suggests strategies to overcome the computational costs due to evaluations of multibody interactions by additional approximation steps with controlled errors. We present numerical examples demonstrating that the inclusion of multibody interactions shows substantial improvement in dynamical simulations, e.g., of rare events and metastability in phase transitions regimes.

Variational Methods for Solving Warped Multirate Partial Differential Algebraic Equations

Roland Pulch

SIAM J. Sci. Comput. 31, pp. 1016-1034 (19 pages)

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
Signals exhibiting amplitude as well as frequency modulation at widely separated time scales arise in radio frequency (RF) applications. A multivariate model yields an adequate representation by decoupling the time scales of involved signals. Consequently, a system of differential algebraic equations (DAEs) modeling the electric circuit changes into a system of partial differential algebraic equations (PDAEs). The determination of an emerging local frequency function is crucial for the efficiency of this approach, since inappropriate choices produce many oscillations in the multivariate solution. Thus, the idea is to reduce oscillating behavior via minimizing the magnitude of partial derivatives. For this purpose, we apply variational calculus to obtain a necessary condition for a specific solution, which represents a minimum of an according functional. This condition can be included in numerical schemes computing the complete solution of the PDAE. Test results confirm that the used strategy ensures an efficient simulation of RF signals.

IDR($s$): A Family of Simple and Fast Algorithms for Solving Large Nonsymmetric Systems of Linear Equations

Peter Sonneveld and Martin B. van Gijzen

SIAM J. Sci. Comput. 31, pp. 1035-1062 (28 pages) | Cited 6 times

Online Publication Date: November 26, 2008

Full Text: | Download PDF

Show Abstract
We present IDR($s$), a new family of efficient, short-recurrence methods for large nonsymmetric systems of linear equations. The new methods are based on the induced dimension reduction (IDR) method proposed by Sonneveld in 1980. IDR($s$) generates residuals that are forced to be in a sequence of nested subspaces. Although IDR($s$) behaves like an iterative method, in exact arithmetic it computes the true solution using at most $N + N/s$ matrix-vector products, with $N$ the problem size and $s$ the codimension of a fixed subspace. We describe the algorithm and the underlying theory and present numerical experiments to illustrate the theoretical properties of the method and its performance for systems arising from different applications. Our experiments show that IDR($s$) is competitive with or superior to most Bi-CG-based methods and outperforms Bi-CGSTAB when $s > 1$.

New Finite Elements for Large-Scale Simulation of Optical Waves

Britta Heubeck, Christoph Pflaum, and Gunther Steinle

SIAM J. Sci. Comput. 31, pp. 1063-1081 (19 pages) | Cited 1 time

Online Publication Date: December 10, 2008

Full Text: | Download PDF

Show Abstract
We present a new method to simulate optical waves in large geometries. This method is based on newly developed finite elements, so-called trigonometric finite wave elements (TFWEs). They are constructed by linear elements as well as by trigonometric functions such that the one-dimensional Helmholtz equation is exactly solved under certain conditions. In comparison with the transfer matrix method, the TFWE method offers equally good results, but it can be extended to higher dimensions and applied to time-dynamic problems. The analysis of TFWEs shows that these elements approximate functions with certain oscillation properties more accurately than standard finite elements. Thus, a finite element discretization with TFWEs leads to a smaller system of equations which eases the solving process. Numerical results obtained by applying the TFWE method to the simulation of the optical wave in distributed feedback lasers are presented.

Augmented Mixed Finite Element Methods for the Stationary Stokes Equations

Leonardo E. Figueroa, Gabriel N. Gatica, and Antonio Márquez

SIAM J. Sci. Comput. 31, pp. 1082-1119 (38 pages) | Cited 1 time

Online Publication Date: December 10, 2008

Full Text: | Download PDF

Show Abstract
In this paper we introduce and analyze two augmented mixed finite element methods for a velocity-pressure-stress formulation of the stationary Stokes equations. Our approach, which extends analog results for linear elasticity problems, is based on the introduction of the Galerkin least-squares-type terms arising from the constitutive and equilibrium equations and the Dirichlet boundary condition for the velocity, all of them multiplied by suitable stabilization parameters. We show that these parameters can be chosen so that the resulting augmented variational formulations are defined by strongly coercive bilinear forms, whence the associated Galerkin schemes become well-posed for any choice of finite element subspaces. In particular, we can use continuous piecewise linear velocities, piecewise constant pressures, and Raviart–Thomas elements for the stresses, thus yielding a number of unknowns behaving asymptotically as 5 times the number of triangles of the triangulation. Alternatively, the above factor reduces to 4 when a second augmented variational formulation, involving only the velocity and the stress as unknowns, is employed. Next, we derive a reliable and efficient residual-based a posteriori error estimator for the augmented mixed finite element schemes. Finally, extensive numerical experiments illustrating the performance of the augmented mixed finite element methods, confirming the properties of the a posteriori estimators, and showing the behavior of the associated adaptive algorithms are reported.

A Three-Dimensional Mixed Finite-Element Approximation of the Semiconductor Energy-Transport Equations

Stephan Gadau and Ansgar Jüngel

SIAM J. Sci. Comput. 31, pp. 1120-1140 (21 pages)

Online Publication Date: December 10, 2008

Full Text: | Download PDF

Show Abstract
The stationary energy-transport equations for semiconductors in three space dimensions are numerically discretized. The physical variables are the electron density, the energy density, and the electric potential. Physically motivated mixed Dirichlet–Neumann boundary conditions are employed. The numerical approximation is based on a hybridized mixed finite-element method with Raviart–Thomas–Nédélec elements, applied to the dual-entropy formulation of the energy-transport model. For the solution of the nonlinear discrete system, a Newton scheme with adaptive potential stepping and two decoupling Gummel-type strategies with reduced rank extrapolation are proposed. Multigate field-effect transistors in two dimensions and three dimensions are numerically simulated.

Asymptotic Stability of a Jump-Diffusion Equation and Its Numerical Approximation

Graeme D. Chalmers and Desmond J. Higham

SIAM J. Sci. Comput. 31, pp. 1141-1155 (15 pages)

Online Publication Date: December 17, 2008

Full Text: | Download PDF

Show Abstract
Asymptotic linear stability is studied for stochastic differential equations (SDEs) that incorporate Poisson-driven jumps and their numerical simulations using theta-method discretizations. The property is shown to have a simple explicit characterization for the SDE, whereas for the discretization a condition is found that is amenable to numerical evaluation. This allows us to evaluate the asymptotic stability behavior of the methods. One surprising observation is that there exist problem parameters for which an explicit, forward Euler method has better stability properties than its trapezoidal and backward Euler counterparts. Other computational experiments indicate that all theta methods reproduce the correct asymptotic linear stability for sufficiently small step sizes. By using a recent result of Appleby, Berkolaiko, and Rodkina, we give a rigorous verification that both linear stability and instability are reproduced for small step sizes. This property is known not to hold for general, nonlinear problems.

Reducing Floating Point Error in Dot Product Using the Superblock Family of Algorithms

Anthony M. Castaldo, R. Clint Whaley, and Anthony T. Chronopoulos

SIAM J. Sci. Comput. 31, pp. 1156-1174 (19 pages)

Online Publication Date: December 17, 2008

Full Text: | Download PDF

Show Abstract
This paper discusses both the theoretical and statistical errors obtained by various well-known dot products, from the canonical to pairwise algorithms, and introduces a new and more general framework that we have named superblock which subsumes them and permits a practitioner to make trade-offs between computational performance, memory usage, and error behavior. We show that algorithms with lower error bounds tend to behave noticeably better in practice. Unlike many such error-reducing algorithms, superblock requires no additional floating point operations and should be implementable with little to no performance loss, making it suitable for use as a performance-critical building block of a linear algebra kernel.

Fast (Structured) Newton Computations

Thomas F. Coleman and Wei Xu

SIAM J. Sci. Comput. 31, pp. 1175-1191 (17 pages)

Online Publication Date: December 17, 2008

Full Text: | Download PDF

Show Abstract
Many vector-valued functions, representing expensive computations, are also structured computations. In this case the calculation of the Newton step can be greatly accelerated by exploiting this structure. It is often not necessary, nor economic, to form the true Jacobian in the process of computing the Newton step; instead, a more cost-effective auxiliary Jacobian matrix is used. This auxiliary matrix can be sparse even when the true Jacobian matrix is dense; consequently, sparse matrix technology can be used, to great speed advantage, both in forming the auxiliary matrix and in solving the auxiliary linear system.

Symmetric Positive Definite Flux-Continuous Full-Tensor Finite-Volume Schemes on Unstructured Cell-Centered Triangular Grids

Helmer André Friis, Michael G. Edwards, and Johannes Mykkeltveit

SIAM J. Sci. Comput. 31, pp. 1192-1220 (29 pages) | Cited 1 time

Online Publication Date: December 31, 2008

Full Text: | Download PDF

Show Abstract
Novel cell-centered full-tensor finite-volume methods are presented for general unstructured grids in two spatial dimensions. The numerical schemes are flux-continuous and based on computing the transmissibilities in a local subcell transform space, ensuring that local flux matrices are symmetric. As a result the global discretization matrix is shown to be symmetric positive definite for any grid type. A symmetric physical space method is also introduced, and the symmetric methods are shown to be closely related. Discrete ellipticity conditions are derived for positive definiteness of the physical space and subcell space schemes. Computational examples are presented for unstructured triangular grids demonstrating good performance of the scheme. The schemes are compared with the so-called multipoint flux approximation (MPFA) O-method [I. Aavatsmark, T. Barkve, Ø. Bøe, and T. Mannseth, SIAM J. Sci. Comput., 19 (1998), pp. 1700–1716]. Good agreement between the methods is obtained, but the new scheme shows improved behavior in challenging cases.

A Hierarchical Multiresolution Adaptive Mesh Refinement for the Solution of Evolution PDEs

S. Jain, P. Tsiotras, and H.-M. Zhou

SIAM J. Sci. Comput. 31, pp. 1221-1248 (28 pages)

Online Publication Date: December 31, 2008

Full Text: | Download PDF

Show Abstract
In this paper, we propose a novel multiresolution adaptive mesh refinement algorithm for solving initial-boundary value problems (IBVP) for evolution PDEs. The proposed algorithm dynamically adapts the grid to any existing or emerging irregularities in the solution, by refining the grid only at those places where the solution exhibits sharp features. The main advantage of the proposed grid adaptation method is that it results in a grid with a fewer number of nodes when compared to adaptive grids generated by existing multiresolution-based mesh refinement techniques. Several examples show the robustness and stability of the proposed algorithm.

Local Discontinuous Galerkin Method for the Hunter–Saxton Equation and Its Zero-Viscosity and Zero-Dispersion Limits

Yan Xu and Chi-Wang Shu

SIAM J. Sci. Comput. 31, pp. 1249-1268 (20 pages)

Online Publication Date: December 31, 2008

Full Text: | Download PDF

Show Abstract
In this paper, we develop, analyze, and test a local discontinuous Galerkin (LDG) method for solving the Hunter–Saxton (HS) equation and its zero-viscosity and zero-dispersion limits. The energy stability for general solutions are proved, and numerical simulation results for different types of solutions of the nonlinear HS equation are provided to illustrate the accuracy and capability of the LDG method. The zero-viscosity and zero-dispersion properties of the HS equation are studied in a numerical simulation.

Accurate Floating-Point Summation Part II: Sign, $K$-Fold Faithful and Rounding to Nearest

Siegfried M. Rump, Takeshi Ogita, and Shin'ichi Oishi

SIAM J. Sci. Comput. 31, pp. 1269-1302 (34 pages) | Cited 3 times

Online Publication Date: December 31, 2008

Full Text: | Download PDF

Show Abstract
In Part II of this paper we first refine the analysis of error-free vector transformations presented in Part I. Based on that we present an algorithm for calculating the rounded-to-nearest result of $s:=\sum p_i$ for a given vector of floating-point numbers $p_i$, as well as algorithms for directed rounding. A special algorithm for computing the sign of $s$ is given, also working for huge dimensions. Assume a floating-point working precision with relative rounding error unit eps. We define and investigate a $K$-fold faithful rounding of a real number $r$. Basically the result is stored in a vector $\mathtt{Res}_{\nu}$ of $K$ nonoverlapping floating-point numbers such that $\sum\mathtt{Res}_{\nu}$ approximates $r$ with relative accuracy $\mathtt{eps}^K$, and replacing $\mathtt{Res}_K$ by its floating-point neighbors in $\sum\mathtt{Res}_{\nu}$ forms a lower and upper bound for $r$. For a given vector of floating-point numbers with exact sum $s$, we present an algorithm for calculating a $K$-fold faithful rounding of $s$ using solely the working precision. Furthermore, an algorithm for calculating a faithfully rounded result of the sum of a vector of huge dimension is presented. Our algorithms are fast in terms of measured computing time because they allow good instruction-level parallelism, they neither require special operations such as access to mantissa or exponent, they contain no branch in the inner loop, nor do they require some extra precision. The only operations used are standard floating-point addition, subtraction, and multiplication in one working precision, for example, double precision. Certain constants used in the algorithms are proved to be optimal.

A Hybrid Phase-Flow Method for Hamiltonian Systems with Discontinuous Hamiltonians

Shi Jin, Hao Wu, and Zhongyi Huang

SIAM J. Sci. Comput. 31, pp. 1303-1321 (19 pages) | Cited 2 times

Online Publication Date: December 31, 2008

Full Text: | Download PDF

Show Abstract
In this paper, we propose a new phase-flow method for Hamiltonian systems with discontinuous Hamiltonians. In the original phase-flow method introduced by Ying and Candès [J. Comput. Phys., 220 (2006), pp. 184–215], the phase map should be smooth to ensure the accuracy of the interpolation. Such an interpolation is inaccurate if the phase map is nonsmooth, for example, when the Hamiltonian is discontinuous. We modify the phase-flow method using a discontinuous Hamiltonian solver and establish the stability (for piecewise constant potentials) of such a solver. This extends the applicability of the highly efficient phase-flow method to singular Hamiltonian systems, with a mild increase of algorithm complexity. Such a particle method can be useful for the computation of high frequency waves through interfaces.

Numerical Integration of Damped Maxwell Equations

M. A. Botchev and J. G. Verwer

SIAM J. Sci. Comput. 31, pp. 1322-1346 (25 pages) | Cited 1 time

Online Publication Date: January 16, 2009

Full Text: | Download PDF

Show Abstract
We study the numerical time integration of Maxwell's equations from electromagnetism. Following the method of lines approach we start from a general semidiscrete Maxwell system for which a number of time-integration methods are considered. These methods have in common an explicit treatment of the curl terms. Central in our investigation is the question how to efficiently raise the temporal convergence order beyond the standard order of two, in particular in the presence of an explicitly or implicitly treated damping term which models conduction.

Existence of Bounded Discrete Steady-State Solutions of the Van Roosbroeck System on Boundary Conforming Delaunay Grids

K. Gärtner

SIAM J. Sci. Comput. 31, pp. 1347-1362 (16 pages) | Cited 1 time

Online Publication Date: January 16, 2009

Full Text: | Download PDF

Show Abstract
The classic van Roosbroeck system describes the carrier transport in semiconductors in a drift diffusion approximation. Its analytic steady state solutions fulfill bounds for some mobility and recombination/generation models. The main goal of this paper is to establish the identical bounds for discrete in space, steady state solutions on 3d boundary conforming Delaunay grids and the classical Scharfetter–Gummel scheme. Together with a uniqueness proof for small applied voltages and the known dissipativity (continuous as well as space and time discrete), these discretization techniques carry over the essential analytic properties to the discrete case. The proofs are of interest for deriving averaging schemes for space or state dependent material parameters, which preserve these qualitative properties, too. To illustrate the properties of the scheme, 1, 4, 16 elementary cells of a modified CoolMOS-like structure are depleted by increasing the applied voltage until steady state avalanche breakdown occurs.

What Makes Molecular Dynamics Work?

Robert D. Skeel

SIAM J. Sci. Comput. 31, pp. 1363-1378 (16 pages) | Cited 3 times

Online Publication Date: January 16, 2009

Full Text: | Download PDF

Show Abstract
The equations of motion for deterministic molecular dynamics (MD) are chaotic, creating problems for their numerical treatment due to the exponential growth of error with time. Indeed, modeling and computational errors overwhelm numerical trajectories in typical simulations. Consequently, accuracy is expected only in a statistical sense, based on random initial conditions. Of great interest then is the relationship between errors in the dynamics and their effects on the accuracy of statistical quantities, specifically, expectations. This article provides a formula for the effect of a perturbation on an ensemble average, which explains the accuracy of such calculations. It also provides a formula for the effect of a perturbation on a time correlation function, which, however, fails to explain accuracy for these calculations. Additionally, this article clarifies the relationships among various dynamical properties of MD and provides an extension to a theory of non-Hamiltonian MD.

Adaptive and Recursive Time Relaxed Monte Carlo Methods for Rarefied Gas Dynamics

Stefano Trazzi, Lorenzo Pareschi, and Bernt Wennberg

SIAM J. Sci. Comput. 31, pp. 1379-1398 (20 pages)

Online Publication Date: January 16, 2009

Full Text: | Download PDF

Show Abstract
Recently a new class of Monte Carlo methods, called time relaxed Monte Carlo (TRMC), designed for the simulation of the Boltzmann equation close to fluid regimes has been introduced [L. Pareschi and G. Russo, SIAM J. Sci. Comput., 23 (2001), pp. 1253–1273]. A generalized Wild sum expansion of the solution is the basis of the simulation schemes. After a splitting of the equation, the time discretization of the collision step is obtained from the Wild sum expansion of the solution by replacing high order terms in the expansion with the equilibrium Maxwellian distribution; in this way speed-up of the methods close to fluid regimes is obtained by efficiently thermalizing particles close to the equilibrium state. In this work we present an improvement of such methods which allows us to obtain an effective uniform accuracy in time without any restriction on the time step and subsequent increase of the computational cost. The main ingredient of the new algorithms is recursivity [L. Pareschi and B. Wennberg, Monte Carlo Methods Appl., 7 (2001), pp. 349–358]. Several techniques can be used to truncate the recursive trees generated by the schemes without degrading the accuracy of the numerical solution. Techniques based on adaptive strategies are presented. Numerical results emphasize the gain of efficiency of the present simulation schemes with respect to standard DSMC (direct simulation Monte Carlo) methods.

A Nonsmooth Multiscale Method for Solving Frictional Two-Body Contact Problems in 2D and 3D with Multigrid Efficiency

Rolf Krause

SIAM J. Sci. Comput. 31, pp. 1399-1423 (25 pages)

Online Publication Date: January 21, 2009

Full Text: | Download PDF

Show Abstract
We present a nonsmooth multiscale method for the numerical solution of frictional contact problems in $2d$ and $3d$. The computational effort is comparable to that of solving linear problems. Our method does not require any regularization, neither of the nonpenetration condition nor of the friction law and can be applied to contact problems with Tresca friction as well as to contact problems with Coulomb friction. For the case of Tresca friction, the global convergence of the method is shown. For the more complicated case of Coulomb friction, we develop a nonsmooth multiscale method which can be directly applied to the corresponding variational quasi-inequality. No outer iteration is required. Moreover, our multiscale approach is general in the sense that it can be used in the context of geometric as well as algebraic multigrid methods. Nonconforming domain decomposition techniques (or mortar) methods are employed in order to enforce the transmission conditions at the interface between different bodies with nonmatching meshes. Numerical examples illustrate the high robustness and efficiency of our method.

Efficient Solvers for a Linear Stochastic Galerkin Mixed Formulation of Diffusion Problems with Random Data

O.G. Ernst, C.E. Powell, D.J. Silvester, and E. Ullmann

SIAM J. Sci. Comput. 31, pp. 1424-1447 (24 pages) | Cited 5 times

Online Publication Date: January 21, 2009

Full Text: | Download PDF

Show Abstract
We introduce a stochastic Galerkin mixed formulation of the steady-state diffusion equation and focus on the efficient iterative solution of the saddle-point systems obtained by combining standard finite element discretizations with two distinct types of stochastic basis functions. So-called mean-based preconditioners, based on fast solvers for scalar diffusion problems, are introduced for use with the minimum residual method. We derive eigenvalue bounds for the preconditioned system matrices and report on the efficiency of the chosen preconditioning schemes with respect to all the discretization parameters.

Constraint Preserving Schemes for Some Gauge Invariant Wave Equations

Snorre H. Christiansen

SIAM J. Sci. Comput. 31, pp. 1448-1469 (22 pages)

Online Publication Date: January 21, 2009

Full Text: | Download PDF

Show Abstract
We study discretizations of the Maxwell–Klein–Gordon equation as an example of a constrained geometric nonlinear evolution partial differential equation. For the temporal gauge we propose a fully discrete scheme which preserves the nonlinear constraint thanks to a special application of Lagrange multipliers. We show that the method generalizes to Hamiltonian wave equations whose kinetic and potential energy are both invariant under a group of transformations, even though the Galerkin spaces are not invariant. We then extend the method to the Lorenz gauge. Numerical results illustrate the discussion.

Coupling of Darcy–Forchheimer and Compressible Navier–Stokes Equations with Heat Transfer

M. Amara, D. Capatina, and L. Lizaik

SIAM J. Sci. Comput. 31, pp. 1470-1499 (30 pages)

Online Publication Date: January 28, 2009

Full Text: | Download PDF

Show Abstract
This paper is devoted to the coupling of a 2D reservoir model with a 1.5D vertical wellbore model, both written in axisymmetric form. The physical problems are, respectively, described by the Darcy–Forchheimer and the compressible Navier–Stokes equations, together with an exhaustive energy equation. Each model was previously studied and its finite element discretization was validated. The two weak problems are bound together by means of transmission conditions at the perforations, yielding a nonstandard mixed formulation. A technical analysis is then carried out and the well-posedness of the time-discretized coupled problem, in both the continuous and the discrete cases, is established. Numerical tests including physical cases are presented, validating the coupled code.

Simulation of Diffraction in Periodic Media with a Coupled Finite Element and Plane Wave Approach

M. Huber, J. Schöberl, A. Sinwel, and S. Zaglmayr

SIAM J. Sci. Comput. 31, pp. 1500-1517 (18 pages) | Cited 1 time

Online Publication Date: January 28, 2009

Full Text: | Download PDF

Show Abstract
The aim of this paper is to discuss simulation methods of diffraction of electromagnetic waves on biperiodic structures. The region with complicated structures is discretized by Nédélec finite elements. In the unbounded homogeneous regions above and below, a plane wave expansion containing the exact far-field pattern is applied. A consistent coupling is achieved by the method of Nitsche. By numerical experiments we investigate the speed of convergence depending on the mesh refinement, the element order, and the number of evanescent waves.

Fourier Analysis of Multigrid Methods on Hexagonal Grids

Guohua Zhou and Scott R. Fulton

SIAM J. Sci. Comput. 31, pp. 1518-1538 (21 pages) | Cited 1 time

Online Publication Date: January 28, 2009

Full Text: | Download PDF

Show Abstract
This paper applies local Fourier analysis to multigrid methods on hexagonal grids. Using oblique coordinates to express the grids and a dual basis for the Fourier modes, the analysis proceeds essentially the same as for rectangular grids. The framework for one- and two-grid analyses is given and then applied to analyze the performance of multigrid methods for the Poisson problem on a hexagonal grid. Numerical results confirm the analysis. Uniform hexagonal grids provide an approximation to spherical geodesic grids; numerical results for the latter show similar performance. While the analysis is similar to that for rectangular grids, the results differ somewhat: full weighting is superior to injection for restriction, Jacobi relaxation performs about as well as Gauss–Seidel relaxation, and underrelaxation is not required for good performance. Also, coarse-fine or four-color ordering (both analogues of red-black ordering on the rectangular grid) improves the performance of Jacobi relaxation, with the latter achieving a smoothing factor of approximately 0.25. An especially simple compact fourth-order discretization works well, and the full multigrid algorithm produces the solution to the level of truncation error in work proportional to the number of unknowns.

Mathematical and Numerical Analysis of a Robust and Efficient Grid Deformation Method in the Finite Element Context

Matthias Grajewski, Michael Köster, and Stefan Turek

SIAM J. Sci. Comput. 31, pp. 1539-1557 (19 pages)

Online Publication Date: January 28, 2009

Full Text: | Download PDF

Show Abstract
Among a variety of grid deformation methods, the method proposed by Bochev, Liao, and de la Pena [Methods Partial Differential Equations, 12 (1996), pp. 489–506], Cai et al. [Comput. Math. Appl., 48 (2004), pp. 1077–1086], and Liao and Anderson [Appl. Anal., 44 (1992), pp. 285–298] is one of the most favorable, because it prevents mesh tangling and offers precise control over the element volumes. Its numerical realization requires only solving a Poisson problem and a system of fully decoupled initial value problems. Many other deformation methods, in contrast, involve the solution of complicated nonlinear partial differential equations (PDEs). In this article, we introduce a generalization of Liao's method, which allows for generating a desired mesh size distribution for quite arbitrary grids without giving rise to mesh tangling. We elaborate on its numerical realization and prove the convergence of our method. Our results are confirmed by numerical experiments.

Solution of One-Time-Step Problems in Elastoplasticity by a Slant Newton Method

P. G. Gruber and J. Valdman

SIAM J. Sci. Comput. 31, pp. 1558-1580 (23 pages)

Online Publication Date: January 28, 2009

Full Text: | Download PDF

Show Abstract
We discuss a solution algorithm for quasi-static elastoplastic problems with hardening. Such problems can be described by a time dependent variational inequality, where the displacement and the plastic strain fields serve as primal variables. After discretization in time, one variational inequality of the second kind is obtained per time step and can be reformulated as a minimization problem with a convex energy functional which depends smoothly on the displacement and nonsmoothly on the plastic strain. There exists an explicit formula for minimizing the energy functional with respect to the plastic strain for a given displacement. By substitution, the energy functional can be written as a functional depending only on the displacement. The theorem of Moreau from convex analysis states that the energy functional is differentiable with an explicitly computable first derivative. The second derivative of the energy functional does not exist, due to the lack of smoothness of the plastic strain across the elastoplastic interface, which separates the continuum in elastically and plastically deformed parts. A Newton-like method exploiting slanting functions of the energy functional's first derivative instead of the nonexistent second derivative is applied. Such a method is called a slant Newton method for short. The local superlinear convergence of the algorithm in the discrete case is shown, and sufficient regularity assumptions are formulated, which would guarantee the local superlinear convergence also in the continuous case.

An Algebraic Treatment of Essential Boundary Conditions in the Particle–Partition of Unity Method

Marc Alexander Schweitzer

SIAM J. Sci. Comput. 31, pp. 1581-1602 (22 pages)

Online Publication Date: February 20, 2009

Full Text: | Download PDF

Show Abstract
This paper is concerned with the treatment of essential boundary conditions in meshfree methods. In particular, we focus on the particle–partition of unity method (PPUM). However, the proposed technique is applicable to any partition of unity-based approach. We present an efficient scheme for the automatic construction of a direct splitting of a PPUM function space into the degrees of freedom suitable for the approximation of the Dirichlet data and the degrees of freedom that remain for the approximation of the PDE by simple linear algebra. Notably, our approach requires no restrictions on the distribution of the discretization points nor on the employed (local) approximation spaces. We attain the splitting of the global function space from the respective direct splittings of the employed local approximation spaces. Hence, the global splitting can be computed with (sub)linear complexity. Due to this direct splitting of the meshfree PPUM function space, we can implement a conforming local treatment of essential boundary data so that the realization of Dirichlet boundary values in the meshfree PPUM is straightforward. The presented approach yields an optimally convergent scheme, which is demonstrated by the presented numerical results.
Close

Your Recent History Recent History Help

Recently Viewed

  1. Iterative Algorithms Based on Decoupling of Deblurring and Denoising for Image Restoration
  2. Constructing Sobol Sequences with Better Two-Dimensional Projections
  3. Asymptotic Sampling Distribution for Polynomial Chaos Representation from Data: A Maximum Entropy an...
  4. Iterated Hard Shrinkage for Minimization Problems with Sparsity Constraints
  5. Continuation of Invariant Subspaces in Large Bifurcation Problems
© 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