SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

BROWSE VOLUMES

Year Range: 

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

Search Issue | RSS Feeds RSS
Next Issue

2012

Volume 34, Issue 1, pp. vii-C41

FREE

SISC Redefined

Hans Petter Langtangen, Editor-in-Chief

SIAM J. Sci. Comput. 34, pp. vii-viii ( pages)

Full Text: | Download PDF

Show Abstract
This issue of SIAM Journal on Scientific Computing marks the beginning of a new era for the journal. During the past decade, the area of scientific computing has continued to evolve and now broadly spans three major areas: the invention and analysis of new core algorithms remain of key importance, but much current research in the field targets challenging scientific or engineering applications, and substantial work goes into advanced software development and exploration of new architectures for large-scale computing.
To properly reflect this ongoing development of the field, the journal is now arranged into three distinct sections to reflect these developments. All three aspects are occasionally involved in today's research projects on advancing computational methods.  It is therefore a natural step for SISC to expand its scope and adapt to properly and accurately reflect changes in the scientific community.
The typical classical SISC paper is now found in the section Methods and Algorithms for Scientific Computing. These papers deal with algorithmic development for scientific or engineering problems of wide interest.  The editorial policy states that “papers in this category may include theoretical analysis, provided that the relevance to applications in science and engineering is demonstrated. They should contain meaningful computational results and theoretical results or strong heuristics supporting the performance of new algorithms."  Accepted papers typically introduce a new algorithm with a performance that is demonstrated to be superior to that of the current state-of-the-art algorithms for the type of problems being addressed.
The new SISC section Computational Methods in Science and Engineering deals with building numerical strategies for exploring challenging scientific or engineering problems through computer simulations. Using the computer as a virtual laboratory in this way is often referred to as computational science and engineering (CSE). Successful papers in this section do not necessarily develop new fundamental algorithms but may offer novelty in the way well-known algorithmic building blocks and computational tools are combined to address a new problem or to explore a previously solved problem in new and superior ways.  To quote the editorial policy: “Papers in this section will typically describe novel methodologies for solving a specific problem in computational science or engineering. They should contain enough information about the application to orient other computational scientists but should omit details of interest mainly to the applications specialist."
The new section on Software and High-Performance Computing (S/HPC) contains papers in which the novelty and innovation occur in the way software is designed and implemented or emerging computing architectures are explored. Significant portions of research budgets in scientific computing projects are spent on these issues, and many projects on scientific discovery through computer models are strongly dependent on innovation in software development, parallel algorithms, and hardware utilization. This section in SISC aims to facilitate the exchange of new ideas and tools and thereby to help raise scientific standards. The editorial policy states: "Papers in this category should concern the development of high-quality computational software, high-performance computing issues, novel architectures, data analysis, or visualization. The primary focus should be on computational methods that have potentially large impact for an important class of scientific or engineering problems."
Authors are asked at submission to choose a SISC section for their manuscript. However, the editors reserve the right to reassign the manuscript to the section they find most appropriate. The editorial structure of SISC has likewise changed to reflect this new structure, and section editors, one for each of the three new sections, have been appointed to help the editor-in-chief implement this new structure.
The new sections in SISC will change neither the purpose nor the quality of the journal. The purpose remains that of advancing computational methods for solving scientific and engineering problems, and the section descriptions above aim to make it clear that SISC papers must address problems and techniques of wide importance, and present new and superior computational approaches. The quality of SISC is widely acknowledged in the scientific community and also reflected in various bibliographic metrics.
SISC dates back to 1980, when the journal was born as SIAM Journal on Scientific and Statistical Computing. Its name was changed to SIAM Journal on Scientific Computing in 1993.  The new section model and the new editorial policy mark yet another important milestone in the journal's history and ensure that SISC will continue to position itself as a leading outlet for research in a rapidly changing and very exciting field.
back to top
RSS Feeds

Wedderburn Rank Reduction and Krylov Subspace Method for Tensor Approximation. Part 1: Tucker Case

S. A. Goreinov, I. V. Oseledets, and D. V. Savostyanov

SIAM J. Sci. Comput. 34, pp. A1-A27 (27 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
New algorithms are proposed for the Tucker approximation of a 3-tensor accessed only through a tensor-by-vector-by-vector multiplication subroutine. In the matrix case, the Krylov methods are methods of choice to approximate the dominant column and row subspaces of a sparse or structured matrix given through a matrix-by-vector operation. Using the Wedderburn rank reduction formula, we propose a matrix approximation algorithm that computes the Krylov subspaces and can be generalized to 3-tensors. The numerical experiments show that on quantum chemistry data the proposed tensor methods outperform the minimal Krylov recursion of Savas and Eldén.

Coupling at a Distance HDG and BEM

Bernardo Cockburn, Francisco-Javier Sayas, and Manuel Solano

SIAM J. Sci. Comput. 34, pp. A28-A47 (20 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
We introduce a new technique for numerically solving exterior Dirichlet boundary-value problems for second-order elliptic equations. It consists of coupling a hybridizable discontinuous Galerkin (HDG) method used for solving the so-called interior problem on a bounded region containing the support of the source term, with a boundary element method (BEM) for solving the problem exterior to that region. The novelty is that the BEM is defined on a suitably chosen, smooth artificial boundary whereas the HDG method is defined on a polyhedral subdomain. Because of the choice of the artificial boundary, we can take advantage of the spectral convergence of the BEM solution and of the simplicity of the corresponding equations. Because the HDG method is defined on a polyhedral subdomain, there is no need to try to fit the mesh to the artificial boundary. Instead, the HDG is coupled at a distance with the BEM by using simple Dirichlet-to-Neumann operators defined in the region between the artificial boundary and the polyhedral subdomain. Numerical experiments displaying the performance of the new technique are presented. Optimal orders of convergence are obtained even though the distance between the artificial boundary and the polyhedral domain is of order $h$.

Modification and Compensation Strategies for Threshold-based Incomplete Factorizations

S. MacLachlan, D. Osei-Kuffuor, and Yousef Saad

SIAM J. Sci. Comput. 34, pp. A48-A75 (28 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
Standard (single-level) incomplete factorization preconditioners are known to successfully accelerate Krylov subspace iterations for many linear systems. The classical modified incomplete LU (MILU) factorization approach improves the acceleration given by (standard) ILU approaches, by modifying the nonunit diagonal in the factorization to match the action of the system matrix on a given vector, typically the constant vector. Here, we examine the role of similar modifications within the dual-threshold ILUT algorithm. We introduce column and row variants of the modified ILUT algorithm and discuss optimal ways of modifying the columns or rows of the computed factors to improve their accuracy and stability. Modifications are considered for both the diagonal and off-diagonal entries of the factors, based on one or many vectors, chosen a priori or through an Arnoldi iteration. Numerical results are presented to support our findings.

New Resolution Strategy for Multiscale Reaction Waves using Time Operator Splitting, Space Adaptive Multiresolution, and Dedicated High Order Implicit/Explicit Time Integrators

Max Duarte, Marc Massot, Stéphane Descombes, Christian Tenaud, Thierry Dumont, Violaine Louvet, and Frédérique Laurent

SIAM J. Sci. Comput. 34, pp. A76-A104 (29 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
We tackle the numerical simulation of reaction-diffusion equations modeling multi-scale reaction waves. This type of problem induces peculiar difficulties and potentially large stiffness which stem from the broad spectrum of temporal scales in the nonlinear chemical source term as well as from the presence of steep spatial gradients in the reaction fronts, spatially very localized. In this paper, we introduce a new resolution strategy based on time operator splitting and space adaptive multiresolution in the context of very localized and stiff reaction fronts. The paper considers a high order implicit time integration of the reaction and an explicit one for the diffusion term in order to build a time operator splitting scheme that exploits efficiently the special features of each problem. Based on recent theoretical studies of numerical analysis such a strategy leads to a splitting time step which is restricted by neither the fastest scales in the source term nor by stability constraints of the diffusive steps but only by the physics of the phenomenon. We aim thus at solving complete models including all time and space scales within a prescribed accuracy, considering large simulation domains with conventional computing resources. The efficiency is evaluated through the numerical simulation of configurations which were so far out of reach of standard methods in the field of nonlinear chemical dynamics for two-dimensional spiral waves and three-dimensional scroll waves, as an illustration. Future extensions of the proposed strategy to more complex configurations involving other physical phenomena as well as optimization capability on new computer architectures are discussed.

Discontinuous Galerkin Approximation of Relaxation Models for Linear and Nonlinear Diffusion Equations

Fausto Cavalli, Giovanni Naldi, and Ilaria Perugia

SIAM J. Sci. Comput. 34, pp. A105-A136 (32 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
In this work we present finite element approximations of relaxed systems for nonlinear diffusion problems, which can also tackle the cases of degenerate and strongly degenerate diffusion equations. Relaxation schemes take advantage of the replacement of the original partial differential equation (PDE) with a semilinear hyperbolic system of equations, with a stiff source term, tuned by a relaxation parameter $\epsilon$. When $\epsilon \rightarrow 0^+$, the system relaxes onto the original PDE: in this way, a consistent discretization of the relaxation system for vanishing $\epsilon$ yields a consistent discretization of the original PDE. The numerical schemes obtained with this procedure do not require solving implicit nonlinear problems and possess the robustness of upwind discretizations. The proposed approximations are based on a discontinuous Galerkin method in space and on suitable implicit-explicit integration in time. Then, in principle, we can achieve any order of accuracy and obtain stable solutions, even when the diffusion equation becomes degenerate and solution singularities develop. Moreover, when needed, we can easily incorporate slope limiters within our schemes in order to handle spurious oscillatory phenomena. Some preliminary theoretical results are given, along with several numerical tests in one and two space dimensions, both for linear and nonlinear diffusion problems, including a degenerate diffusion equation, that provide numerical evidence of the properties of the presented approach.

The Direct Richardson $p$th Order (DRp) Schemes: A New Class of Time Integration Schemes for Stochastic Differential Equations

Pavel P. Popov and Stephen B. Pope

SIAM J. Sci. Comput. 34, pp. A137-A160 (24 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
We describe a new family of weak $p$th order accurate SDE time integration schemes, called the direct Richardson $p$th order accurate (DRp) schemes. The DRp schemes use the idea of Richardson extrapolation on Euler time steps, performed by way of an acceptance-rejection algorithm. Previous applications of Richardson extrapolation to the Euler scheme are applicable only when the objective is to estimate a functional of the final distribution of the process. In contrast, provided that the diffusion matrix is strictly positive definite, the DRp class of schemes can be used in all applications which require a weak SDE time integration scheme. Numerical results have been obtained, and a comparison is made between the second- and third-order accurate DRp schemes and other modern SDE time integration schemes, indicating that the DRp schemes incur less error than standard algorithms based on Ito–Taylor expansions, and have similar computational efficiency. Finally, we provide a proof of the convergence properties of the DRp schemes.

State Trajectory Compression for Optimal Control with Parabolic PDEs

Martin Weiser and Sebastian Götschel

SIAM J. Sci. Comput. 34, pp. A161-A184 (24 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
In optimal control problems with nonlinear time-dependent three-dimensional (3D) PDEs, full four-dimensional (4D) discretizations are usually prohibitive due to the storage requirement. For this reason gradient- and Newton-type methods working on the reduced functional are often employed. The computation of the reduced gradient requires one solve of the state equation forward in time and one backward solve of the adjoint equation. The state enters into the adjoint equation, again requiring the storage of a full 4D data set. We propose to use lossy compression algorithms, using an inexact but cheap predictor for the state data, with additional entropy coding of prediction errors. As the data is used inside a discretized, iterative algorithm, lossy compression maintaining a certain error bound turns out to be sufficient. We provide an error analysis, derive expected compression rates, and present numerical examples validating the results.

Adaptive and Stochastic Algorithms for Electrical Impedance Tomography and DC Resistivity Problems with Piecewise Constant Solutions and Many Measurements

Kees van den Doel and Uri M. Ascher

SIAM J. Sci. Comput. 34, pp. A185-A205 (21 pages) | Cited 1 time

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
This article develops fast numerical methods for the practical solution of the famous electrical impedance tomography and DC resistivity problems in the presence of discontinuities and potentially many experiments or data. Based on a Gauss–Newton (GN) approach coupled with preconditioned conjugate gradient (PCG) iterations, we propose two algorithms. One determines adaptively the number of inner PCG iterations required to stably and effectively carry out each GN iteration. The other algorithm, useful especially in the presence of many experiments, employs a randomly chosen subset of experiments at each GN iteration that is controlled using a cross validation approach. Numerical examples demonstrate the efficacy of our algorithms.

Communication-optimal Parallel and Sequential QR and LU Factorizations

James Demmel, Laura Grigori, Mark Hoemmen, and Julien Langou

SIAM J. Sci. Comput. 34, pp. A206-A239 (34 pages)

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
We present parallel and sequential dense QR factorization algorithms that are both optimal (up to polylogarithmic factors) in the amount of communication they perform and just as stable as Householder QR. We prove optimality by deriving new lower bounds for the number of multiplications done by “non-Strassen-like” QR, and using these in known communication lower bounds that are proportional to the number of multiplications. We not only show that our QR algorithms attain these lower bounds (up to polylogarithmic factors), but that existing LAPACK and ScaLAPACK algorithms perform asymptotically more communication. We derive analogous communication lower bounds for LU factorization and point out recent LU algorithms in the literature that attain at least some of these lower bounds. The sequential and parallel QR algorithms for tall and skinny matrices lead to significant speedups in practice over some of the existing algorithms, including LAPACK and ScaLAPACK, for example, up to 6.7 times over ScaLAPACK. A performance model for the parallel algorithm for general rectangular matrices predicts significant speedups over ScaLAPACK.

A Matrix-free Approach for Solving the Parametric Gaussian Process Maximum Likelihood Problem

Mihai Anitescu, Jie Chen, and Lei Wang

SIAM J. Sci. Comput. 34, pp. A240-A262 (23 pages)

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
Gaussian processes are the cornerstone of statistical analysis in many application areas. Nevertheless, most of the applications are limited by their need to use the Cholesky factorization in the computation of the likelihood. In this work, we present a matrix-free approach for computing the solution of the maximum likelihood problem involving Gaussian processes. The approach is based on a stochastic programming reformulation followed by sample average approximation applied to either the maximization problem or its optimality conditions. We provide statistical estimates of the approximate solution. The method is illustrated on several examples where the data is provided on a regular or irregular grid. In the latter case, the action of a covariance matrix on a vector is computed by means of fast multipole methods. For each of the examples presented, we demonstrate that the approach scales linearly with an increase in the number of sites.

Compact and Stable Discontinuous Galerkin Methods for Convection-Diffusion Problems

S. Brdar, A. Dedner, and R. Klöfkorn

SIAM J. Sci. Comput. 34, pp. A263-A282 (20 pages)

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
We present a new scheme, the compact discontinuous Galerkin 2 (CDG2) method, for solving nonlinear convection-diffusion problems together with a detailed comparison to other well-accepted DG methods. The new CDG2 method is similar to the CDG method that was recently introduced in the work of Perraire and Persson for elliptic problems. One main feature of the CDG2 method is the compactness of the stencil which includes only neighboring elements, even for higher order approximation. Theoretical results showing coercivity and stability of CDG2 and CDG for the Poisson and the heat equation are given, providing computable bounds on any free parameters in the scheme. In numerical tests for an elliptic problem, a scalar convection-diffusion equation, and for the compressible Navier–Stokes equations, we demonstrate that the CDG2 method slightly outperforms similar methods in terms of $L^2$-accuracy and CPU time.

Enhancing Quasi-Monte Carlo Methods by Exploiting Additive Approximation for Problems in Finance

Xiaoqun Wang

SIAM J. Sci. Comput. 34, pp. A283-A308 (26 pages)

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
The valuation of many financial securities can be formulated as high-dimensional integrals for which quasi-Monte Carlo (QMC) methods are becoming indispensable numerical tools. It is known that some typical high-dimensional problems in finance and in many other fields have a strong additive nature. This paper investigates the enhancement of QMC methods by exploiting the additive approximation in combination with suitable transformation methods. To this end, we study (a) the approximation of a multivariate function by a sum of one-dimensional functions; and (b) the effective use of the possible strong additive property inherent in finance problems to enhance QMC. For problem (a), we establish a relationship between the approximation quality resulting from anchored decomposition and the global sensitivity indices. For problem (b), we propose new methods for efficiency improvement—the additive weighted uniform sampling and the additive control variate, which exactly exploit the strong additive structure. The required importance distribution and control variate are found constructively based on the anchored additive approximation with a suitably chosen anchor. Numerical examples on bond valuation and option pricing illustrate that the proposed methods reduce the variance by huge factors in QMC. Transformation methods play a crucial role, since they may enhance the degree of additivity of a function. Traditional methods, such as Brownian bridge and principal component analysis, do not necessarily work well, whereas methods that take into account the underlying functions are robust and powerful.

Wavelet Collocation Method and Multilevel Augmentation Method for Hammerstein Equations

Hideaki Kaneko, Khomsan Neamprem, and Boriboon Novaprateep

SIAM J. Sci. Comput. 34, pp. A309-A338 (30 pages)

Online Publication Date: February 07, 2012

Full Text: | Download PDF

Show Abstract
A wavelet collocation method for nonlinear Hammerstein equations is formulated. A sparsity in the Jacobian matrix is obtained which gives rise to a fast algorithm for nonlinear solvers such as the Newton's method and the quasi-Newton method. A fast multilevel augmentation method is developed on a transformed nonlinear equation which gives an additional saving in computational time.

NR$xx$ Simulation of Microflows with Shakhov Model

Zhenning Cai, Ruo Li, and Zhonghua Qiao

SIAM J. Sci. Comput. 34, pp. A339-A369 (31 pages)

Online Publication Date: February 07, 2012

Full Text: | Download PDF

Show Abstract
In this paper, we propose a method to simulate the microflows with Shakhov model using the NR$xx$ method developed in [Z. Cai and R. Li, SIAM J. Sci. Comput., 32 (2010), pp. 2875–2907; Z. Cai, R. Li, and Y. Wang, Commun. Comput. Phys., 11 (2012), pp. 1415–1438; Z. Cai, R. Li, and Y. Wang, J. Sci. Comput., to appear]. The equation under consideration is the Boltzmann equation with force terms, and the Shakhov model is adopted to achieve the correct Prandtl number. As the focus of this paper, we derive a uniform framework for different order moment systems on the wall boundary conditions, which is a major difficulty in the moment methods. Numerical examples for both steady and unsteady problems are presented to show the convergence in the number of moments.

Time Implicit High-Order Discontinuous Galerkin Method with Reduced Evaluation Cost

Florent Renac, Claude Marmignon, and Frédéric Coquel

SIAM J. Sci. Comput. 34, pp. A370-A394 (25 pages)

Online Publication Date: February 07, 2012

Full Text: | Download PDF

Show Abstract
An efficient and robust time integration procedure for a high-order discontinuous Galerkin method is introduced for solving nonlinear second-order partial differential equations. The time discretization is based on an explicit formulation for the hyperbolic term and an implicit formulation for the parabolic term. The procedure uses an iterative algorithm with reduced evaluation cost. The size of the linear system to be solved is greatly reduced thanks to partial uncoupling in space between low-order and high-order degrees of freedom. Numerical examples are presented for the nonlinear convection-diffusion equation in one and two dimensions including steady and unsteady flow problems. The performance of the present method is investigated in terms of CPU time and compared to a fully implicit method. A von Neumann stability analysis is carried out in order to determine the stability and damping properties of the method. Besides a fairly reduced CPU effort, numerical results demonstrate better convergence properties of the present algorithm when compared to the fully implicit method.

On Finding Multiple Solutions to a Singularly Perturbed Neumann Problem

Ziqing Xie, Yongjun Yuan, and Jianxin Zhou

SIAM J. Sci. Comput. 34, pp. A395-A420 (26 pages)

Online Publication Date: February 07, 2012

Full Text: | Download PDF

Show Abstract
In this paper, in order to numerically solve for multiple positive solutions to a singularly perturbed Neumann boundary value problem in mathematical biology and other applications, a local minimax method is modified with new local mesh refinement and other strategies. Algorithm convergence and other related properties are verified. Motivated by the numerical algorithm and convinced by the numerical results, a Morse index approach is used to identify the Morse index of the root solution $u^1_\varepsilon=1$ at any perturbation value, its bifurcation points, and then the critical perturbation value. Many interesting numerical solutions are computed for the first time and displayed with their contours and mesh profiles to illustrate the theory and method.

The Optimized Schwarz Method with a Coarse Grid Correction

Olivier Dubois, Martin J. Gander, Sébastien Loisel, Amik St-Cyr, and Daniel B. Szyld

SIAM J. Sci. Comput. 34, pp. A421-A458 (38 pages)

Online Publication Date: February 09, 2012

Full Text: | Download PDF

Show Abstract
Optimized Schwarz methods (OSMs) use Robin transmission conditions across the subdomain interfaces. The Robin parameter can then be optimized to obtain the fastest convergence. A new formulation is presented with a coarse grid correction. The optimal parameter is computed for a model problem on a cylinder, together with the corresponding convergence factor which is smaller than that of classical Schwarz methods. A new coarse space is presented, suitable for OSM. Numerical experiments illustrating the effectiveness of OSM with a coarse grid correction, both as an iteration and as a preconditioner, are reported.

Deflation, Projector Preconditioning, and Balancing in Iterative Substructuring Methods: Connections and New Results

Axel Klawonn and Oliver Rheinbach

SIAM J. Sci. Comput. 34, pp. A459-A484 (26 pages)

Online Publication Date: February 09, 2012

Full Text: | Download PDF

Show Abstract
In this paper, projector preconditioning, also known as the deflation method, as well as the balancing preconditioner are applied to the dual-primal finite element tearing and interconnecting (FETI-DP) and balancing domain decomposition by constraints (BDDC) methods in order to create a second, independent coarse problem. This may help to extend the parallel scalability of classical FETI-DP and BDDC methods without the use of inexact solvers and may also be used to improve the robustness, e.g., for almost incompressible elasticity problems. Connections of FETI-DP methods applying a transformation of basis using a larger coarse space with a corresponding FETI-DP method using projector preconditioning or balancing are pointed out. It is then shown that the methods have essentially the same spectrum. Numerical results for compressible and almost incompressible linear elasticity are provided. The sensitivity of the projection methods to an inexact computation of the projections is numerically investigated and a different behavior for projector preconditioning and the balancing preconditioner is found.

Computing All or Some Eigenvalues of Symmetric $\mathcal{H}_{\ell}$-Matrices

Peter Benner and Thomas Mach

SIAM J. Sci. Comput. 34, pp. A485-A496 (12 pages)

Online Publication Date: February 14, 2012

Full Text: | Download PDF

Show Abstract
We use a bisection method [B. Parlett, The Symmetric Eigenvalue Problem, Prentice-Hall, Englewood Cliffs, NJ, 1980, p. 51] to compute the eigenvalues of a symmetric $\mathcal{H}_{\ell}$-matrix $M$. The number of negative eigenvalues of $M-\mu I$ is computed via the LDL$^T$ factorization of $M-\mu I$. For dense matrices, the LDL$^T$ factorization is too expensive to yield an efficient eigenvalue algorithm in general, but not so for $\mathcal{H}_{\ell}$-matrices. In the special structure of $\mathcal{H}_{\ell}$-matrices there is an LDL$^T$ factorization with linear-polylogarithmic complexity. The bisection method requires only matrix-size independent many iterations to find an eigenvalue up to the desired accuracy, so that an eigenvalue can be found in linear-polylogarithmic time. For all $n$ eigenvalues, $\mathcal{O}(n^{2}(\log n )^{4}\log( \Vert {M}\Vert_{2}/\epsilon_{\text{ev}}))$ flops are needed to compute all eigenvalues with an accuracy $\epsilon_{\text{ev}}$. It is also possible to compute only eigenvalues in a specific interval or the $j$th smallest one. Numerical experiments demonstrate the efficiency of the algorithm, in particular for the case when some interior eigenvalues are required. (A corrected PDF is attached to this article.)

Solving Dirichlet Boundary-value Problems on Curved Domains by Extensions from Subdomains

Bernardo Cockburn and Manuel Solano

SIAM J. Sci. Comput. 34, pp. A497-A519 (23 pages)

Online Publication Date: February 14, 2012

Full Text: | Download PDF

Show Abstract
We present a technique for numerically solving Dirichlet boundary-value problems for second-order elliptic equations on domains $\Omega$ with curved boundaries. This is achieved by using suitably defined extensions from polyhedral subdomains $\textsf{D}_h$; the problem of dealing with curved boundaries is thus reduced to the evaluations of simple line integrals. The technique is independent of the representation of the boundary, of the space dimension, and allows for the use of only polyhedral elements. We apply this technique to the hybridizable discontinuous Galerkin method and provide extensive numerical experiments showing that the convergence properties of the resulting method are the same as those for the case in which $\Omega=\textsf{D}_h$ whenever the distance of $\textsf{D}_h$ to $\partial\Omega$ is of order $h$. In particular, we find that when using polynomial approximations of degree $k$ on the polyhedral subdomains $\textsf{D}_h$, both the scalar and vector unknowns converge with the optimal order of $k+1$ in the whole domain $\Omega$ for $k\ge0$. Moreover, we also find that for $k\ge1$, both a postprocessing of the scalar unknown in the polyhedral subdomain $\textsf{D}_h$ as well as its extension from it converge with order $k+2$. The extension also converges with order $2$ for $k=0$.

Maxwellian Decay for Well-balanced Approximations of a Super-characteristic Chemotaxis Model

Laurent Gosse

SIAM J. Sci. Comput. 34, pp. A520-A545 (26 pages) | Cited 1 time

Online Publication Date: February 23, 2012

Full Text: | Download PDF

Show Abstract
We focus on the numerical simulation of a one-dimensional model of chemotaxis dynamics (proposed by Greenberg and Alt [Trans. Amer. Math. Soc., 300 (1987), pp. 235–258] in a bounded domain by means of a previously introduced well-balanced (WB) and asymptotic-preserving (AP) scheme [L. Gosse, J. Math. Anal. Appl., (2011)]. We are especially interested in studying the decay onto numerical steady-states for two reasons: (1) conventional upwind schemes have been shown to stabilize onto spurious non-Maxwellian states (with a very big mass flux; see, e.g., [F. Guarguaglini et al., Discrete Contin. Dyn. Syst. Ser. B, 12 (2009), pp. 39–76]) and (2) the initial data lead to a dynamic which is mostly super-characteristic in the sense of [S. Jin and M. Katsoulakis, SIAM J. Appl. Math., 61 (2000), pp. 273–292]; thus the stability results of Gosse do not apply. A reflecting boundary condition which is compatible with the well-balanced character is presented; a mass-preservation property is proved and some results on super-characteristic relaxation are recalled. Numerical experiments with coarse computational grids are presented in detail: they illustrate the bifurcation diagrams of Guarguaglini et al., which relate the total initial mass of cells with the time-asymptotic values of the chemoattractant substance on each side of the domain. It is shown that the WB scheme stabilizes correctly onto zero-mass flow rate (hence Maxwellian) steady-states which agree with the aforementioned bifurcation diagrams. The evolution in time of residues is commented for every considered test case.
back to top
RSS Feeds

Controlling Errors in Recursive Fermi–Dirac Operator Expansions with Applications in Electronic Structure Theory

Emanuel H. Rubensson

SIAM J. Sci. Comput. 34, pp. B1-B23 (23 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
Recursive Fermi–Dirac operator expansion is an efficient way to compute one-electron density matrices in electronic structure theory. The convergence is rapid and depends only weakly on the conditioning of the problem and, for many systems, the computational cost increases only linearly with system size. In this article, errors introduced when evaluating the recursive expansion are analyzed and schemes to control the forward error are proposed. The error has previously been analyzed for explicit schemes working at zero electronic temperature [J. Chem. Phys., 128 (2008), 074106]. Here, implicit schemes [Phys. Rev. B, 68 (2003), 233104] working at zero or finite temperature are treated. The proposed schemes for error control are demonstrated by tight-binding as well as density functional theory electronic structure calculations on several test systems. Condition numbers for the problem of computing the density matrix are derived, giving quantitative insight into under what circumstances a temperature dependent formulation results in better conditioning. It is shown that for the considered recursive expansions, the number of matrix-matrix multiplications needed to compute the density matrix increases only with the squared logarithm of the condition number of the problem.

Computing Stochastic Traveling Waves

G. J. Lord and V. Thümmler

SIAM J. Sci. Comput. 34, pp. B24-B43 (20 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
The aim of this paper is to investigate new numerical methods for computing traveling wave solutions and new ways for estimating characteristic properties such as wave speed for stochastically forced partial differential equations. As a particular example we consider the Nagumo equation with multiplicative noise which we mainly consider in the Stratonovich sense. A standard approach for determining the position and hence speed of a wave is to compute the evolution of a level set. We compare this approach against an alternative where the wave position is found by minimizing the $L^2$ norm against a fixed reference profile. This approach can be used to freeze (or stop) the wave and obtain a stochastic partial differential algebraic equation that we then discretize and solve. Although attractive because it leads to a smaller domain size, it can be numerically unstable due to large convection terms. We compare numerically the different approaches for estimating the wave speed. We use these techniques to investigate the effect of both Itô and Stratonovich noise on the Nagumo equation as correlation length and noise intensity increase.

Uncertainty Quantification given Discontinuous Model Response and a Limited Number of Model Runs

Khachik Sargsyan, Cosmin Safta, Bert Debusschere, and Habib Najm

SIAM J. Sci. Comput. 34, pp. B44-B64 (21 pages)

Online Publication Date: February 02, 2012

Full Text: | Download PDF

Show Abstract
We outline a methodology for forward uncertainty quantification in systems with uncertain parameters, discontinuous model response, and a limited number of model runs. Our approach involves two stages. First we detect the discontinuity with Bayesian inference, thus obtaining a probabilistic representation of the discontinuity curve for arbitrarily distributed input parameters. Then, employing the Rosenblatt transform, we construct spectral representations of the uncertain model output, using polynomial chaos (PC) expansions on either side of the discontinuity curve, leading to an averaged PC representation of the forward model response that allows efficient uncertainty quantification. We obtain PC modes by either orthogonal projection or Bayesian inference, and argue for a hybrid approach that targets a balance between the accuracy provided by the orthogonal projection and the flexibility provided by the Bayesian inference. The uncertain model output is then computed by taking an ensemble average over PC expansions corresponding to sampled realizations of the discontinuity curve. The methodology is demonstrated on synthetic examples of discontinuous model response with adjustable sharpness and structure.

Computation of Tail Probabilities via Extrapolation Methods and Connection with Rational and Padé Approximants

Philippe Gaudreau, Richard M. Slevinsky, and Hassan Safouhi

SIAM J. Sci. Comput. 34, pp. B65-B85 (21 pages)

Online Publication Date: February 09, 2012

Full Text: | Download PDF

Show Abstract
We use the recently developed algorithm for the $G_{n}^{(1)}$ transformation to approximate tail probabilities of the normal distribution, the gamma distribution, the student's $t$-distribution, the inverse Gaussian distribution, and Fisher's $F$ distribution. Using this algorithm, which can be computed recursively when using symbolic programming languages, we are able to compute these integrals to high predetermined accuracies. Previous to this contribution, the evaluation of these tail probabilities using the $G_{n}^{(1)}$ transformation required symbolic computation of large determinants. With the use of our algorithm, the $G_n^{(1)}$ transformation can be performed relatively easily to produce explicit approximations. After a brief theoretical study, a connection between the $G_n^{(1)}$ transformation and rational and Padé approximants is established.

A Generalized Mixed Hybrid Mortar Method for Solving Flow in Stochastic Discrete Fracture Networks

G. Pichot, J. Erhel, and J.-R. de Dreuzy

SIAM J. Sci. Comput. 34, pp. B86-B105 (20 pages)

Online Publication Date: February 28, 2012

Full Text: | Download PDF

Show Abstract
The simulation of flow in fractured media requires handling both a large number of fractures and a complex interconnecting network of these fractures. Networks considered in this paper are three-dimensional domains made up of two-dimensional fractures intersecting each other and randomly generated. Due to the stochastic generation of fractures, intersections can be highly intricate. The numerical method must generate a mesh and define a discrete problem for any discrete fracture network (DFN). A first approach [Erhel, de Dreuzy, and Poirriez, SIAM J. Sci. Comput., 31 (2009), pp. 2688–2705] is to generate a conforming mesh and to apply a mixed hybrid finite element method. However, the resulting linear system becomes very large when the network contains many fractures. Hence a second approach [Pichot, Erhel, and de Dreuzy, Appl. Anal., 89 (2010), pp. 1629–1643] is to generate a nonconforming mesh, using an independent mesh generation for each fracture. Then a Mortar technique applied to the mixed hybrid finite element method deals with the nonmatching grids. When intersections do not cross or overlap, pairwise Mortar relations for each intersection are efficient [Pichot, Erhel, and de Dreuzy, 2010]. But for most random networks, discretized intersections involve more than two fractures. In this paper, we design a new method generalizing the previous one and that is applicable for stochastic networks. The main idea is to combine pairwise Mortar relations with additional relations for the overlapping part. This method still ensures the continuity of fluxes and heads and still yields a symmetric positive definite linear system. Numerical experiments show the efficiency of the method applied to complex stochastic fracture networks. We also study numerical convergence when reducing the mesh step. This method makes it easy to perform mesh optimization and appears to be a very promising tool to simulate flow in multiscale fracture networks.
back to top
RSS Feeds

Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques

Hasan Metin Aktulga, Sagar A. Pandit, Adri C. T. van Duin, and Ananth Y. Grama

SIAM J. Sci. Comput. 34, pp. C1-C23 (23 pages)

Online Publication Date: January 31, 2012

Full Text: | Download PDF

Show Abstract
Modeling atomic and molecular systems requires computation-intensive quantum mechanical methods such as, but not limited to, density functional theory [R. A. Friesner, Proc. Natl. Acad. Sci. USA, 102 (2005), pp. 6648–6653]. These methods have been successful in predicting various properties of chemical systems at atomistic scales. Due to the inherent nonlocality of quantum mechanics, the scalability of these methods ranges from O($N^3$) to O($N^7$) depending on the method used and approximations involved. This significantly limits the size of simulated systems to a few thousand atoms, even on large scale parallel platforms. On the other hand, classical approximations of quantum systems, although computationally (relatively) easy to implement, yield simpler models that lack essential chemical properties such as reactivity and charge transfer. The recent work of van Duin et al. [J. Phys. Chem. A, 105 (2001), pp. 9396–9409] overcomes the limitations of nonreactive classical molecular dynamics (MD) approximations by carefully incorporating limited nonlocality (to mimic quantum behavior) through an empirical bond order potential. This reactive classical MD method, called ReaxFF, achieves essential quantum properties, while retaining the computational simplicity of classical MD, to a large extent. Implementation of reactive force fields presents significant algorithmic challenges. Since these methods model bond breaking and formation, efficient implementations must rely on complex dynamic data structures. Charge transfer in these methods is accomplished by minimizing electrostatic energy through charge equilibration. This requires the solution of large linear systems ($10^8$ degrees of freedom and beyond) with shielded electrostatic kernels at each time-step. Individual time-steps are themselves typically in the range of tenths of femtoseconds, requiring optimizations within and across time-steps to scale simulations to nanoseconds and beyond, where interesting phenomena may be observed. In this paper, we present implementation details of sPuReMD (serial Purdue reactive molecular dynamics program), a unique reactive classical MD code. We describe various data structures, and the charge equilibration solver at the core of the simulation engine. This Krylov subspace solver relies on a preconditioner based on incomplete LU factorization with thresholds (ILUT), specially targeted to our application. We comprehensively validate the performance and accuracy of sPuReMD on a variety of hydrocarbon systems. In particular, we show excellent per-time-step time, linear time scaling in system size, and a low memory footprint. sPuReMD is a freely distributed software with GPL and is currently being used to model diverse systems ranging from oxidative stress in biomembranes to strain relaxation in Si-Ge nanorods.

Framework for Massively Parallel Adaptive Finite Element Computational Fluid Dynamics on Tetrahedral Meshes

Niclas Jansson, Johan Hoffman, and Johan Jansson

SIAM J. Sci. Comput. 34, pp. C24-C41 (18 pages)

Online Publication Date: February 16, 2012

Full Text: | Download PDF

Show Abstract
In this paper we describe a general adaptive finite element framework for unstructured tetrahedral meshes without hanging nodes suitable for large scale parallel computations. Our framework is designed to scale linearly to several thousands of processors, using fully distributed and efficient algorithms. The key components of our implementation, local mesh refinement and load balancing algorithms, are described in detail. Finally, we present a theoretical and experimental performance study of our framework, used in a large scale computational fluid dynamics computation, and we compare scaling and complexity of different algorithms on different massively parallel architectures.
Close

Your Recent History Recent History Help

Recently Viewed

  1. A Bootstrap Algebraic Multilevel Method for Markov Chains
  2. Simulation, Characterization, and Optimization of Metabolic Models with the High Performance Systems...
  3. A General Framework for Deriving Integral Preserving Numerical Methods for PDEs
  4. Tensor-Train Decomposition
  5. On Moduli of Rings and Quadrilaterals: Algorithms and Experiments
© 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