SIAM Digital Library
 
 
 

SIAM J. on Scientific Computing

Year Range: 

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

Search Issue | RSS Feeds RSS
Previous Issue

2008

Volume 30, Issue 6, pp. vii-3288

* Special Issue on Computational Science and Engineering

back to top
RSS Feeds

Special Issue on Computational Science and Engineering

Chris Johnson, David Keyes, and Ulrich Ruede

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

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Leading-edge science and engineering depend on advanced computing for understanding, prediction, and control. In response to these needs, the field of computational science and engineering (CS&E) is evolving rapidly, to the point that it is now widely considered to be a new discipline by itself and a third pillar of the scientific enterprise, a peer alongside theory and physical experiment. CS&E is unique in that it enables progress in virtually all other disciplines by providing windows of discovery when traditional means of research reach their limits.
Because of its flexibility, computer simulation has become a universal tool. A simulation may serve as a virtual microscope that lets scientists observe the world of quantum physics much smaller than an atom, or it may be employed as a virtual telescope that allows us to explore how galaxies are forming in the universe. In this way, CS&E helps scientists to reach beyond our physical limitations in space and time. When physical experiments are too dangerous, too time-consuming, too expensive, or simply impossible, then advanced simulation techniques enable us to perform virtual experiments and to obtain scientific results otherwise beyond our capabilities.
Just as scientists employ computer simulations to improve our understanding of the world, engineers can use them to design better products. Virtual prototypes—simulated devices and processes—can increasingly replace their real counterparts to save time and to reduce cost, while simulation-based optimization techniques help engineers to develop better products. CS&E is used for research into innovative materials, to make vehicles safer and environmentally friendlier, to construct buildings that can withstand earthquakes and hurricanes, to improve energy efficiency everywhere, and to advance health care, among many other applications.
To begin CS&E research into a specific problem, researchers generally build a simulator that can correctly represent what has been experimentally observed. The creation of a successful simulator is often a very useful accomplishment in itself, since the work required for such a creation will typically lead to a much improved understanding of the physics (or chemistry, etc.) behind the observations. A working simulation will illuminate which effects are essential in the physics under study, which are essential for the observed behavior, and which can be neglected. In this phase, relatively simple scenarios and model geometries will often be sufficient to perform the necessary work. For example, it may be that a simple two-dimensional model will enable researchers to explore the essential mechanism under consideration. Up to this level, qualitative results are often sufficient.
However, as computational modeling continues to develop, CS&E research will aim to develop computational models that are sophisticated enough to be predictive. Good predictions will require quantitative accuracy, as well as the ability to address more complicated scenarios. This is often very difficult in practice. For example, methods that have been proven useful for computing the flow in rectangular driven cavities may still be very far from suitable in predicting the flow in an aneurysm in the human brain. Data structures and boundary conditions along with data input and output become much more challenging at this stage.
Once CS&E research has progressed to the level of being predictive, the final and ultimate step it must take to become indispensable in an engineering application is optimization. In general, an engineer wants to understand and make predictions about a device because he or she wants to improve the design; thus, he or she will want to use computing to compare different design choices and to decide how an engineering problem may be solved in the best possible way. Too often, practice in CS&E is still that a human varies the parameters, and predictive simulation is used to evaluate how well design criteria have been met. CS&E must enable the researcher to perform these optimization tasks seamlessly within a single model.
In addition to application domain knowledge, CS&E research requires expertise in advanced computing (elements of applied mathematics, scientific computing, and computer science). Additionally, many CS&E tasks involve complex data analysis to provide suitable input for the computations and visualizations required to help humans interpret the simulation outcome. Thus, many CS&E problems can be characterized by a pipeline that includes modeling techniques (mathematical and geometric), simulation techniques (discretizations, algorithms, data structures, software frameworks, and problem solving environments), and analysis techniques (data mining, data management, visualization, as well as error, sensitivity, stability, and uncertainty analyses).
This special issue on CS&E of the SIAM Journal on Scientific Computing is the result of a special call for papers that was announced at the SIAM Conference on Computational Science and Engineering held in Costa Mesa, California, February 19–23, 2007. Eighty-nine papers were submitted to the special issue and subjected to a review process as rigorous as the standard review procedures of the SIAM Journal on Scientific Computing. Three guest editors-in-chief headed an editorial board of 22 well-known researchers in CS&E. After a peer review process with two or more anonymous referees for each paper, 26 papers were eventually selected for publication and are now collected in this volume. The editors would herewith like to thank all contributors, including those whose papers do not appear in this volume. Several of these submissions may eventually appear as regular papers in SISC or other journals. As a result of our strict selectivity, we believe that this volume may be considered a milestone of case studies marking the current state of the art for research in CS&E by example.
As these papers clearly illustrate, CS&E research does not primarily seek to achieve new results for a given domain science but rather to provide new computational tools to researchers. Even though these tools are regularly developed to address challenging problems in specific domain sciences, they are generally transferable to other applications. For example, a method to simulate large ensembles of interacting particles may be useful to an astrophysicist studying galaxy formation as well as to a nanotechnology researcher exploring molecular dynamics. If the methods were not transferable to other applications, a collection of papers of such wide span of different applications into a single volume would be of dubious value. Ultimately, CS&E aims at developing a universal set of simulation methods and tools for the scientists and engineers of the future.
Nevertheless, high-level CS&E research differs from research in mathematics or computer science in that it is oriented at one or more realistic applications of a domain science or engineering discipline. Such problems often involve complicated three-dimensional geometries, multiple interacting scales, heterogeneities, anisotropies, and multiphysical or biological descriptions. Thus, they often thwart rigorous proofs of accuracy or efficiency, and so CS&E research must regularly address validation and verification by other means than traditional mathematical analysis.
Additionally, CS&E research cannot be limited to any single step of the simulation pipeline. A mathematical simulation algorithm is useful in CS&E only when it is also implemented on a real computer, and a new parallel processing paradigm becomes relevant for CS&E only when it is demonstrated to be applicable for a large-scale numerical computation. Of course, research for an individual piece of the CS&E pipeline, such as the development of a new finite element method, may also be very useful for progress in CS&E, but taken by itself such a development remains merely work in numerical analysis. In other words, CS&E is not just a new name for existing computational disciplines, nor—least of all—is it just mathematics in which the proofs are half finished. In fact, as opposed to other mathematical sciences, CS&E often achieves its progress through a clever combination of techniques and methods employed for the different stages of the CS&E pipeline. In such a case, the innovation of excellent research may consist in the creativity needed to synthesize a computational solution for a complex problem from the right building blocks.
Since this is a genuine volume on CS&E, it is impossible to classify the papers in it into single subtopics of CS&E. Just as genuine CS&E research must synthesize its solution from various building blocks, so each of these papers addresses several aspects of what we consider to be typical of CS&E research. Therefore, the assignment of a paper to any single such category is necessarily somewhat arbitrary. However, all papers of this special issue share a focus on the development of innovative computational methodology in the context of challenging applications, in the spirit of our rationale above.
Flow problems, when interpreted in a wide sense, form the largest class of applications that occurs in this special issue. This section includes papers on iterative solvers for flow problems, such as “Newton-GMRES Preconditioning for Discontinuous Galerkin Discretizations of the Navier–Stokes Equations” by Persson and Peraire and turbulence modeling in “A Galerkin-Characteristic Method for Large-Eddy Simulation of Turbulent Flow and Heat Transfer” by El Amrani and Seaid. The flow applications in this volume span the scales, from climate modeling in “A New Approach to Atmospheric General Circulation Model: Global Cloud Resolving Model NICAM and Its Computational Performance” by Tomita, Goto, and Satoh to microfluidics in “Numerical Simulation of Particle Transport in a Drift Ratchet” by Brenk, Bungartz, Mehl, Muntean, Neckel, and Weinzierl. Also dealing with flow problems are “Hybrid Simulations of Reaction-Diffusion Systems in Porous Media” by Tartakovsky, Tartakovsky, Scheibe, and Meakin and “Output Functional Control for Nonlinear Equations Driven by Anisotropic Mesh Adaptation: The Navier–Stokes Equations” by Micheletti and Perotto. Reflecting a current trend in computational methods research, four papers in this volume deal with biomedical applications. They include “A Parallel Multilevel Technique for Solving the Bidomain Equation on a Human Heart with Purkinje Fibers and a Torso Model” by Washio, Okada, and Hisada; “Bridging Scales: A Three-Dimensional Electromechanical Finite Element Model of Skeletal Muscle” by Roehrle, Davidson, and Pullan; “Parallel Minimum $p$-Norm Solution of the Biomagnetic Inverse Problem for Realistic Signals Using Exact Hessian-Vector Products” by Bücker, Buecker, and Rupp; and “Long-Time Simulations on High Resolution Meshes to Model Calcium Waves in a Heart Cell” by Gobbert.
Adaptive methods in which the evolving solution and some aspect of the method (typically the discretization, but increasingly algorithmic parameters as well) interact so that the latter may be optimized, have been and clearly will remain a hot topic in CS&E research. Papers such as “PML Enhanced with a Self-Adaptive Goal-Oriented $hp$-Finite Element Method: Simulation of Through-Casing Borehole Resistivity Measurements” by Pardo, Demkowicz, Torres-Verdín, and Michler; “A Framework for the Adaptive Finite Element Solution of Large-Scale Inverse Problems” by Bangerth; and “Adaptive Discrete Galerkin Methods Applied to the Chemical Master Equation” by Deuflhard, Huisinga, Jahnke, and Wulkow illustrate such adaptive techniques in complex applications.
Reflecting a growth area whose techniques now have a home in a specialized SIAM journal, several papers are concerned with image processing, including “Adaptive Mesh Refinement for Nonparametric Image Registration” by Haber, Heldmann, and Modersitzki; “A Variational Shape Optimization Approach for Image Segmentation with the Mumford–Shah Functional” by Dogan, Morin, and Nochetto; and “Brain–Tumor Interaction Biophysical Models for Medical Image Registration” by Hogea, Davatzikos, and Biros.
Two papers deal with fast transforms, including “An Adaptive Multilevel Wavelet Solver for Elliptic Equations on an Optimal Spherical Geodesic Grid” by Mehra and Kevlahan and “Fast Three-Dimensional Discrete Cosine Transform” by Lee, Chan, and Adjeroh. The paper by Narumi, Kameoka, Taiji, and Yasuoka presents the innovative use of multicore processors for molecular dynamics under the title “Accelerating Molecular Dynamics Simulations on PlayStation 3 Platform Using Virtual-GRAPE Programming Model.” The focus is on stochastic modeling and simulation methods in two papers: “Fast Monte Carlo Simulation Methods for Biological Reaction-Diffusion Systems in Solution and on Surfaces” by Kerr, Bartol, Kaminsky, Dittrich, Chang, Baden, Sejnowski, and Stiles; and “Towards a Statistical Theory of Texture Evolution in Polycrystals” by Barmak, Emelianenko, Golovaty, Kinderlehrer, and Ta'asan.
Finally, a group of five papers can be counted under the topic of optimization, identification and inverse problems (though these topics of course also play a role in the imaging papers above and others). The papers here include “Numerical Solution of an Inverse Problem of Imaging of Antipersonnel Land Mines by the Globally Convergent Convexification Algorithm” by Xin and Klibanov; “ORBIT: Optimization by Radial Basis Function Interpolation in Trust-Regions” by Wild, Regis, and Shoemaker; “Multigrid One-Shot Method for State Constrained Aerodynamic Shape Optimization” by Hazra; “Incremental Identification of Transport Coefficients in Convection-Diffusion Systems” by Karalashvili, Gross, Mhamdi, Reusken, and Marquardt; and “Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space” by Bui-Thanh, Willcox, and Ghattas.
To conclude, the guest editors-in-chief would like to thank again all contributors, associate editors, and the anonymous referees from industry, government laboratories, and universities across the globe for the work and enthusiasm they have put into this publication project. As SIAM seeks to chart a path on the research frontier of CS&E that is responsive to the interests of its members and illuminating to a vast community of developers and users not yet SIAM members, this collection helps to define CS&E, to transmit exciting developments, and above all to improve and focus our collective conversation about CS&E. What is the future role of CS&E in SIAM publications? This special issue complements a growing book series in leading to answers!

Newton-GMRES Preconditioning for Discontinuous Galerkin Discretizations of the Navier–Stokes Equations

P.-O. Persson and J. Peraire

SIAM J. Sci. Comput. 30, pp. 2709-2733 (25 pages) | Cited 2 times

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
We study preconditioners for the iterative solution of the linear systems arising in the implicit time integration of the compressible Navier–Stokes equations. The spatial discretization is carried out using a discontinuous Galerkin method with fourth order polynomial interpolations on triangular elements. The time integration is based on backward difference formulas resulting in a nonlinear system of equations which is solved at each timestep. This is accomplished using Newton's method. The resulting linear systems are solved using a preconditioned GMRES iterative algorithm. We consider several existing preconditioners such as block Jacobi and Gauss–Seidel combined with multilevel schemes which have been developed and tested for specific applications. While our results are consistent with the claims reported, we find that these preconditioners lack robustness when used in more challenging situations involving low Mach numbers, stretched grids, or high Reynolds number turbulent flows. We propose a preconditioner based on a coarse scale correction with postsmoothing based on a block incomplete LU factorization with zero fill-in (ILU0) of the Jacobian matrix. The performance of the ILU0 smoother is found to depend critically on the element numbering. We propose a numbering strategy based on minimizing the discarded fill-in in a greedy fashion. The coarse scale correction scheme is found to be important for diffusion dominated problems, whereas the ILU0 preconditioner with the proposed ordering is effective at handling the convection dominated case. While little can be said in the way of theoretical results, the proposed preconditioner is shown to perform remarkably well for a broad range of representative test problems. These include compressible flows ranging from very low Reynolds numbers to fully turbulent flows using the Reynolds averaged Navier–Stokes equations discretized on highly stretched grids. For low Mach number flows, the proposed preconditioner is more than one order of magnitude more efficient than the other preconditioners considered.

A Galerkin-Characteristic Method for Large-Eddy Simulation of Turbulent Flow and Heat Transfer

Mofdi El-Amrani and Mohammed Seaid

SIAM J. Sci. Comput. 30, pp. 2734-2754 (21 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
This work aims at the development of a nonoscillatory Galerkin-characteristic method for large-eddy simulation of turbulent flow and heat transfer. The method is based on combining the modified method of characteristics with a Galerkin finite element discretization of the incompressible Navier–Stokes/Boussinesq equations in primitive variables. It can be interpreted as a fractional step technique where the convective part and the Stokes/Boussinesq part are treated separately. A limiting procedure is implemented for the reconstruction of numerical solutions at the departure points. The main feature of the proposed Galerkin-characteristic method is that, due to the Lagrangian treatment of convection, the standard Courant–Friedrichs–Levy condition is relaxed, and the time truncation errors are reduced in the Stokes/Boussinesq part. To solve the generalized Stokes/Boussinesq problem we implement a conjugate gradient algorithm. This method avoids projection techniques and does not require any special correction for the pressure. We verify the method for an advection-diffusion equation with a known analytical solution and for the benchmark problem of mixed convection flow in a squared cavity. We also present numerical results for a problem of heat transport in the Strait of Gibraltar. The Galerkin-characteristic method has been found to be feasible and satisfactory.

A New Approach to Atmospheric General Circulation Model: Global Cloud Resolving Model NICAM and its Computational Performance

Hirofumi Tomita, Koji Goto, and Masaki Satoh

SIAM J. Sci. Comput. 30, pp. 2755-2776 (22 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
We have developed a new type of atmospheric general circulation model for super–high-resolution climate studies. To avoid the pole problem, the model uses a modified icosahedral grid and a newly developed nonhydrostatic scheme. In this paper, we concentrate on the computational aspects of the model, including optimization for a massively parallel supercomputer. We optimized the code to achieve fast calculation on the Earth Simulator (ES) and compared the computational performance with that of the spectral transform model. Our results show that our new approach may be a prototype for future atmospheric general circulation models.

Numerical Simulation of Particle Transport in a Drift Ratchet

Markus Brenk, Hans-Joachim Bungartz, Miriam Mehl, Ioan L. Muntean, Tobias Neckel, and Tobias Weinzierl

SIAM J. Sci. Comput. 30, pp. 2777-2798 (22 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
The directed transport of microparticles depending on their size is the basis for particle sorting methods that are of utmost importance in, for example, life sciences. A drift ratchet is a Brownian motor that allows for such a directed transport. Hereby, the particle motion is induced by a combination of the Brownian motion and asymmetries stemming, for example, from the domain's geometry, electrical fields, or transient pressure boundary conditions. We simulate a particular drift ratchet which consists of a matrix of pores with asymmetrically oscillating diameter wherein a fluid with suspended particles is pumped forward and backward, and where the particles' long-term transport direction depends on their size. Thus, this setup allows for continuous and parallel particle separation, which has been shown experimentally. However, for a deeper understanding and for an optimized parameters' choice, further investigations, i.e., simulations, are necessary. In this paper, we present first results achieved with our parallel three-dimensional simulation codes applied to a still simplified scenario. This simplification is necessary to isolate different phenomena (e.g., asymmetries and Brownian motion) to check their relevance for the particle transport. The simulation codes are based on (adaptive) Cartesian grids in combination with finite volume and finite element discretizations. Cartesian grids allow for a very efficient implementation of the solver algorithms and an efficient balanced parallelization via domain decomposition. The achieved simulation results show the effectiveness of our approach and give some strong hints on a directed particle transport already with the simplified model we used here.

Hybrid Simulations of Reaction-Diffusion Systems in Porous Media

A. M. Tartakovsky, D. M. Tartakovsky, T. D. Scheibe, and P. Meakin

SIAM J. Sci. Comput. 30, pp. 2799-2816 (18 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Hybrid or multiphysics algorithms provide an efficient computational tool for combining micro- and macroscale descriptions of physical phenomena. Their use becomes imperative when microscale descriptions are too computationally expensive to be conducted in the whole domain, while macroscale descriptions fail in a small portion of the computation domain. We present a hybrid algorithm to model a general class of reaction-diffusion processes in granular porous media, which includes mixing-induced mineral precipitation on, or dissolution of, the porous matrix. These processes cannot be accurately described using continuum (Darcy-scale) models. The pore-scale/Darcy-scale hybrid is constructed by coupling solutions of the reaction-diffusion equations (RDE) at the pore scale with continuum Darcy-level solutions of the averaged RDEs. The resulting hybrid formulation is solved numerically by employing a multiresolution meshless discretization based on the smoothed particle hydrodynamics method. This ensures seamless noniterative coupling of the two components of the hybrid model. Computational examples illustrate the accuracy and efficiency of the hybrid algorithm.

Output Functional Control for Nonlinear Equations Driven by Anisotropic Mesh Adaption: The Navier–Stokes Equations

Stefano Micheletti and Simona Perotto

SIAM J. Sci. Comput. 30, pp. 2817-2854 (38 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
The contribution of this paper is twofold. Firstly, moving from the very well-known dual-weighted residual (DWR) method, we set up a theoretical framework for a goal-oriented a posteriori analysis of nonlinear partial differential equations accounting for different approximations of the primal and dual problems; nonhomogeneous Dirichlet boundary conditions, even different on passing from the primal to the dual problem; the error due to data approximation; and the effect of a possible stabilization. Secondly, moving from this framework and employing anisotropic interpolation error estimates, a sound anisotropic mesh adaption procedure is devised for the numerical approximation of the Navier–Stokes equations by continuous piecewise linear finite elements. The resulting adaptive procedure is thoroughly addressed and validated on some relevant test cases.

A Parallel Multilevel Technique for Solving the Bidomain Equation on a Human Heart with Purkinje Fibers and a Torso Model

Takumi Washio, Jun-ichi Okada, and Toshiaki Hisada

SIAM J. Sci. Comput. 30, pp. 2855-2881 (27 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
In this paper, we present a multigrid method and its implementation on parallel computers to solve the bidomain equation that appears in excitation propagation analysis of the human heart with the torso. The bidomain equation is discretized with the finite element method on a composite mesh composed of a fine voxel mesh around the heart and a coarse voxel mesh covering the torso. The extracellular potential problem on the torso is formulated as a variational problem with a constraint at the interface of the fine and coarse meshes. We show that this formulation naturally satisfies the conservation property of the electric currents and fits into the multilevel adaptive solution technique framework. We also present our special treatment of the Purkinje fiber network in the multigrid algorithm where it is modeled as multiway branching lines connected to the nodes in the voxel mesh of the heart. A parallel implementation of the proposed multigrid algorithm on distributed memory computers is presented and the performance is evaluated using real-life applications.

Bridging Scales: A Three-Dimensional Electromechanical Finite Element Model of Skeletal Muscle

Oliver Röhrle, John B. Davidson, and Andrew J. Pullan

SIAM J. Sci. Comput. 30, pp. 2882-2904 (23 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
This paper introduces a framework for skeletal muscles that couples outputs from a detailed biophysically based electrophysiological cell model to a three-dimensional continuum-based finite element model of muscle mechanics. Due to the unique manner in which a skeletal muscle is activated, specifically the fact that neighboring fibers are electrically isolated and can act independently of each other, a completely new and novel coupling framework has been created. Within this framework, the electrical activity within a fiber is modeled with a biophysically based cell model, which is itself an amalgamation of several existing cell models. From this amalgamated cell model, specific output parameters that describe the level of crossbridge activity are computed and stored within a lookup table. This lookup table is then used to map the appropriate level of activity to all fibers within the muscle. To link the level of activity to a three-dimensional finite flement model of a skeletal muscle, which is based on principles of continuum mechanics, an upscaling method is introduced to compensate for the fact that the finite element mesh does not attempt to separately represent each individual fiber. This upscaling method allows the stress equilibrium equations to be computed at each Gauss point based on different values of the cell model outputs in all the neighboring cells. Since adjacent fibers can operate independently, the cell model outputs used in the finite element solution of the finite elasticity equations are discontinuous. The behavior and performance of the entire coupling framework is carefully analyzed in some simple test cases analyzing the reduction of the discretization error with respect to a sequence of uniformly refined meshes and different activation patterns. The results show that the error-reduction factors obtained from the electromechanical framework using triquadratic Lagrange and tricubic Hermite basis functions in solving the Galerkin finite element stress equilibrium equations are very similar to those obtained from a mechanics-only continuum-based model. Following this, an example of this process applied to the lateral pterygoid muscle is presented. The proposed framework can be used, for example, to investigate the mechanical effects with respect to cellular changes or to analyze the effects of different neuromuscular activation patterns on the tissue response.

Parallel Minimum $p$-Norm Solution of the Neuromagnetic Inverse Problem for Realistic Signals Using Exact Hessian-Vector Products

H. Martin Bücker, Roland Beucker, and André Rupp

SIAM J. Sci. Comput. 30, pp. 2905-2921 (17 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
In the neuromagnetic inverse problem, one is interested in determining the current density inside the human brain from measurements of the magnetic field recorded outside the head. From a numerical point of view, the solution of this inverse problem is challenging not only in terms of nonuniqueness and time complexity but also with respect to numerical stability. An efficient and robust computational technique is presented that finds the minimum $p$-norm solution of the neuromagnetic inverse problem. The approach is based on carefully combining a subspace trust-region algorithm for the solution of an unconstrained nonlinear optimization problem, automatic differentiation for the truncation-error free evaluation of first- and second-order derivatives, and shared-memory parallelization using the OpenMP programming paradigm. Using actual measurements obtained from a head phantom model as well as realistic data sets of middle-latency auditory evoked field data, it is demonstrated that a valid reconstruction of neuromagnetic activity is achieved for values of $p$ less than 2.

Long-Time Simulations on High Resolution Meshes to Model Calcium Waves in a Heart Cell

Matthias K. Gobbert

SIAM J. Sci. Comput. 30, pp. 2922-2947 (26 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
A model for the flow of calcium on the scale of one heart cell is given by a system of time-dependent reaction-diffusion equations coupled by nonlinear reaction terms. Calcium ions enter into the cell at release units distributed throughout the cell and then diffuse. At each release unit, the probability for calcium to be released increases along with the concentration of calcium, thus creating a feedback loop of waves regenerating themselves repeatedly. The validation of this model requires simulations on the time scale of several repeated waves and on the spatial scale of the entire cell. This requires long-time studies on spatial meshes that need to have a high resolution to resolve the positions of the calcium release units throughout the entire cell. We detail the development of a special-purpose numerical method and parallel implementation for this problem. Parallel performance studies demonstrate the scalability of the implementation on a distributed-memory cluster with low-latency interconnect. Convergence studies verify convergence to analytical expectations and confirm the appropriateness of all numerical parameters. Application studies on the desired time and length scales confirm that the model exhibits the desired feedback mechanism for calcium currents through the release units at suitable high levels, but the long-time studies demonstrate also that the current model with its present parameters leads to excessive calcium concentrations over time. This phenomenon could only be observed using a computational method able to reach laboratory scale final times for a domain on the scale of a complete cell.

PML Enhanced with a Self-Adaptive Goal-Oriented $hp$-Finite Element Method: Simulation of Through-Casing Borehole Resistivity Measurements

D. Pardo, L. Demkowicz, C. Torres-Verdín, and C. Michler

SIAM J. Sci. Comput. 30, pp. 2948-2964 (17 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
We describe the application of a perfectly matched layer (PML) combined with a self-adaptive goal-oriented $hp$-finite element (FE) method to the simulation of borehole resistivity measurements. The adaptive refinements and fast convergence of the self-adaptive $hp$-FE method enhance the performance of the PML, thereby enabling the accurate and efficient truncation of the computational domain in open-domain problems. We apply this method to the simulation of axisymmetric through-casing resistivity measurements acquired in a borehole environment that are typically used for the assessment of rock formation properties. Our numerical results confirm the accuracy and efficiency of our method and provide evidence of highly accurate and reliable simulations of borehole logging measurements in the presence of a conductive steel casing and material contrast of fourteen orders of magnitude in electrical conductivity. Moreover, the combination of adaptivity and PML enables us to significantly reduce the size of the computational domain with substantial savings in computer time and memory.

A Framework for the Adaptive Finite Element Solution of Large-Scale Inverse Problems

Wolfgang Bangerth

SIAM J. Sci. Comput. 30, pp. 2965-2989 (25 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Since problems involving the estimation of distributed coefficients in partial differential equations are numerically very challenging, efficient methods are indispensable. In this paper, we will introduce a framework for the efficient solution of such problems. This comprises the use of adaptive finite element schemes, solvers for the large linear systems arising from discretization, and methods to treat additional information in the form of inequality constraints on the parameter to be recovered. The methods to be developed will be based on an all-at-once approach, in which the inverse problem is solved through a Lagrangian formulation. The main feature of the paper is the use of a continuous (function space) setting to formulate algorithms, in order to allow for discretizations that are adaptively refined as nonlinear iterations proceed. This entails that steps such as the description of a Newton step or a line search are first formulated on continuous functions and only then evaluated for discrete functions. On the other hand, this approach avoids the dependence of finite dimensional norms on the mesh size, making individual steps of the algorithm comparable even if they used differently refined meshes. Numerical examples will demonstrate the applicability and efficiency of the method for problems with several million unknowns and more than 10,000 parameters.

Adaptive Discrete Galerkin Methods Applied to the Chemical Master Equation

P. Deuflhard, W. Huisinga, T. Jahnke, and M. Wulkow

SIAM J. Sci. Comput. 30, pp. 2990-3011 (22 pages) | Cited 2 times

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
In systems biology, the stochastic description of biochemical reaction kinetics is increasingly being employed to model gene regulatory networks and signaling pathways. Mathematically speaking, such models require the numerical solution of the underlying evolution equation, known as the chemical master equation (CME). Until now, the CME has primarily been treated by Monte Carlo techniques, the most prominent of which is the stochastic simulation algorithm [D. T. Gillespie, J. Comput. Phys., 22 (1976), pp. 403–434]. The paper presents an alternative, which focuses on the discrete partial differential equation (PDE) structure of the CME. This allows us to adopt ideas from adaptive discrete Galerkin methods as first suggested by Deuflhard and Wulkow [IMPACT Comput. Sci. Engrg., 1 (1989), pp. 269–301] for polyreaction kinetics and independently developed by Engblom. From the two different options for discretizing the CME as a discrete PDE, Engblom chose the method of lines approach (first space, then time), whereas we strongly advocate use of the Rothe method (first time, then space) for clear theoretical and algorithmic reasons. Numerical findings at two rather challenging problems illustrate the promising features of the proposed method and, at the same time, indicate lines of necessary further improvement of the method worked out here.

Adaptive Mesh Refinement for Nonparametric Image Registration

Eldad Haber, Stefan Heldmann, and Jan Modersitzki

SIAM J. Sci. Comput. 30, pp. 3012-3027 (16 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Three-dimensional (3D) image registration is a computationally intensive problem which is commonly solved in medical imaging. The complexity of the problem stems from its size and nonlinearity. In this paper we present an approach that drastically reduces the problem size by using adaptive mesh refinement. Our approach requires special and careful discretization of the variational form on adaptive quad/octree grids. It further requires an appropriate refinement criterion. We show that in some cases this approach can reduce the computational time by a factor of approximately 10 or so in two dimensions and 5 in three dimensions compared to the nonadaptive approach.

A Variational Shape Optimization Approach for Image Segmentation with a Mumford–Shah Functional

Günay Doğan, Pedro Morin, and Ricardo H. Nochetto

SIAM J. Sci. Comput. 30, pp. 3028-3049 (22 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
We introduce a novel computational method for a Mumford–Shah functional, which decomposes a given image into smooth regions separated by closed curves. Casting this as a shape optimization problem, we develop a gradient descent approach at the continuous level that yields nonlinear PDE flows. We propose time discretizations that linearize the problem and space discretization by continuous piecewise linear finite elements. The method incorporates topological changes, such as splitting and merging for detection of multiple objects, space–time adaptivity, and a coarse-to-fine approach to process large images efficiently. We present several simulations that illustrate the performance of the method and investigate the model sensitivity to various parameters.

Brain–Tumor Interaction Biophysical Models for Medical Image Registration

Cosmina Hogea, Christos Davatzikos, and George Biros

SIAM J. Sci. Comput. 30, pp. 3050-3072 (23 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
State-of-the art algorithms for deformable image registration are based on the minimization of an image similarity functional that is regularized by adding a penalty term on the deformation map. The penalty function typically represents a smoothness regularization. In this article, we use a constrained optimization formulation in which the image similarity functional is coupled to a biophysical model. This formulation is pertinent when the data have been generated by imaging tissue that undergoes deformations due to an actual biophysical phenomenon. Such is the case of coregistering tumor-bearing brain images from the same individual. We present an approximate model that couples tumor growth with the mechanical deformations of the surrounding brain tissue. We consider primary brain tumors—in particular, gliomas. Glioma growth is modeled by a reaction-advection-diffusion PDE, with a two-way coupling with the underlying tissue elastic deformation. Tumor bulk, infiltration, and subsequent mass effects are not regarded separately but are captured by the model itself in the course of its evolution. Our formulation allows for updating the tumor diffusion coefficient following structural displacements caused by tumor growth/infiltration. Our forward problem implementation builds on the PETSc library of Argonne National Laboratory. Our reformulation results in a very small parameter space, and we use the derivative-free optimization library APPSPACK of Sandia National Laboratories. We test the forward model and the optimization framework by using landmark-based similarity functions and by applying it to brain tumor data from clinical and animal studies. State-of-the-art registration algorithms fail in such problems due to excessive deformations. We compare our results with previous work in our group, and we observed up to 50% improvement in landmark deformation prediction. We present preliminary validation results in which we were able to reconstruct deformation fields using four degrees of freedom. Our study demonstrates the validity of our formulation and points to the need for richer datasets and fast optimization algorithms.

An Adaptive Multilevel Wavelet Solver for Elliptic Equations on an Optimal Spherical Geodesic Grid

Mani Mehra and Nicholas K.-R. Kevlahan

SIAM J. Sci. Comput. 30, pp. 3073-3086 (14 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
An adaptive multilevel wavelet solver for elliptic equations on an optimal spherical geodesic grid is developed. The method is based on second-generation spherical wavelets on almost uniform optimal spherical geodesic grids. It is an extension of the adaptive multilevel wavelet solver [O. V. Vasilyev and N. K.-R. Kevlahan, J. Comput. Phys., 206 (2005), pp. 412–431] to curved manifolds. Wavelet decomposition is used for grid adaption and interpolation. A hierarchical finite difference scheme based on the wavelet multilevel decomposition is used to approximate the Laplace–Beltrami operator. The optimal spherical geodesic grid [Internat. J. Comput. Geom. Appl., 16 (2006), pp. 75–93] is convergent in terms of local mean curvature and has lower truncation error than conventional spherical geodesic grids. The overall computational complexity of the solver is $O(\mathcal{N})$, where $\mathcal{N}$ is the number of grid points after adaptivity. The accuracy and efficiency of the method is demonstrated for the spherical Poisson equation. Although the present paper considers the sphere, the strength of this new method is that it can be extended easily to other curved manifolds by choosing an appropriate coarse approximation and using recursive surface subdivision.

Fast Three-Dimensional Discrete Cosine Transform

M. C. Lee, Raymond K. W. Chan, and Don A. Adjeroh

SIAM J. Sci. Comput. 30, pp. 3087-3107 (21 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
A fast three-dimensional discrete cosine transform algorithm (3D FCT) and a fast 3D inverse cosine transform (3D IFCT) algorithm are presented, suitable for analysis of 3D data points. Many existing algorithms for three-dimensional data points make use of either the 1D cosine transform or both the 2D and 1D cosine transforms. Existing algorithms based on the 1D discrete cosine transform (DCT) apply the separable 1D transform to the data points in the $x$, $y$, and $z$ directions, respectively, while those based on 2D and 1D transforms apply the 2D cosine transform for the $x$-$y$ planes and then the 1D cosine transform in the $z$ direction. The proposed 3D DCT algorithms handle the 3D data points directly and have been shown to be computationally efficient. They involve a 3D decomposition process where a data volume is recursively decomposed in each dimension until unit data cubes are obtained. The algorithms are presented in the form of a signal flow graph which captures the various computations involved. A complexity analysis along with empirical results is included, demonstrating the performance of the proposed direct 3D DCT algorithms. As 3D FCT and IFCT are symmetric and relatively fast, they can be used in any application requiring a real-time symmetric codec, such as video conferencing, online multiparty video games, and three-dimensional graphics rendering.

Accelerating Molecular Dynamics Simulations on PlayStation 3 Platform Using Virtual-GRAPE Programming Model

Tetsu Narumi, Shun Kameoka, Makoto Taiji, and Kenji Yasuoka

SIAM J. Sci. Comput. 30, pp. 3108-3125 (18 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Molecular dynamics (MD) simulation requires huge computational power because each atom interacts with another by long-range forces, such as Coulomb and van der Waals forces. Therefore, parallel computers or special accelerators, such as MDGRAPE-3, are required for accelerating MD simulations. A video game processor in a Sony PlayStation 3 or NVIDIA's graphic accelerator card is also a candidate hardware for accelerating MD simulations, because the peak performance of the latest video game processors exceeds that of a current PC's CPU, and they are also very cost-effective. However, the software development for these processors requires much more time compared to CPUs of PCs because the hardware has a lot of parallel processing. We propose the virtual-GRAPE programming model to utilize the hardware resource of video game processors with minimum time for software development. GRAPE is a special-purpose computer used to accelerate particle-based simulations: astrophysical or MD simulations. Under the Virtual-GRAPE model, the subroutine whose calculation speed is accelerated by the special hardware, GRAPE, is replaced with a specially tuned subroutine to be used without the accelerator. We implemented this model in a PlayStation 3 to accelerate the “sander” MD module in the AMBER software package. We were able to achieve an acceleration of 20 times, compared to a serial job using an Intel Xeon 5160 processor. Its performance cost is far superior to that of a PC or an MDGRAPE-3. To obtain the highest performance from the subroutine, most of the arithmetic operations in the tuned routine were performed with single precision accuracy, which is sufficient for MD simulations.

Fast Monte Carlo Simulation Methods for Biological Reaction-Diffusion Systems in Solution and on Surfaces

Rex A. Kerr, Thomas M. Bartol, Boris Kaminsky, Markus Dittrich, Jen-Chien Jack Chang, Scott B. Baden, Terrence J. Sejnowski, and Joel R. Stiles

SIAM J. Sci. Comput. 30, pp. 3126-3149 (24 pages) | Cited 2 times

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Many important physiological processes operate at time and space scales far beyond those accessible to atom-realistic simulations, and yet discrete stochastic rather than continuum methods may best represent finite numbers of molecules interacting in complex cellular spaces. We describe and validate new tools and algorithms developed for a new version of the MCell simulation program (MCell3), which supports generalized Monte Carlo modeling of diffusion and chemical reaction in solution, on surfaces representing membranes, and combinations thereof. A new syntax for describing the spatial directionality of surface reactions is introduced, along with optimizations and algorithms that can substantially reduce computational costs (e.g., event scheduling, variable time and space steps). Examples for simple reactions in simple spaces are validated by comparison to analytic solutions. Thus we show how spatially realistic Monte Carlo simulations of biological systems can be far more cost-effective than often is assumed, and provide a level of accuracy and insight beyond that of continuum methods.

Towards a Statistical Theory of Texture Evolution in Polycrystals

K. Barmak, M. Emelianenko, D. Golovaty, D. Kinderlehrer, and S. Ta'asan

SIAM J. Sci. Comput. 30, pp. 3150-3169 (20 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
Most technologically useful materials possess polycrystalline microstructures composed of a large number of small monocrystalline grains separated by grain boundaries. The energetics and connectivity of the grain boundary network play a crucial role in determining the properties of a material across a wide range of scales. A central problem in materials science is to develop technologies capable of producing an arrangement of grains—a texture—that provides for a desired set of material properties. One of the most challenging aspects of this problem is to understand the role of topological reconfigurations during coarsening. Here we propose an upscaling procedure suitable for large complex systems. The procedure is based on numerical experimentation combined with stochastic tools and consists of large-scale numerical simulations of a system at a microscopic level, statistical analysis of the microscopic data, and formulation of the model based on stochastic characteristics predicted by the statistical analysis. This approach promises to be valuable in establishing the effective model of microstructural evolution in realistic two- and three-dimensional systems. To test ideas we use our upscaling procedure to study the mesoscopic behavior of a reduced one-dimensional network of grain boundaries. Despite the simplicity of its formulation, this model exhibits highly nontrivial behavior characterized by growth and disappearance of grain boundaries and develops probability distributions similar to those observed in higher-dimensional simulations. Here we focus on the grain deletion events which are common to all coarsening systems.

Numerical Solution of an Inverse Problem of Imaging of Antipersonnel Land Mines By the Globally Convergent Convexification Algorithm

Jianguo Xin and Michael V. Klibanov

SIAM J. Sci. Comput. 30, pp. 3170-3196 (27 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
The globally convergent convexification algorithm has been successfully applied to some coefficient inverse problems. In this paper a new numerical implementation of this algorithm is presented and applied to a model problem of imaging of antipersonnel plastic land mines using electric signal propagation. A smoothing technique is introduced to smooth out the spikes on the image for the reconstructed unknown coefficient. Some artifacts of resulting images are explained based on numerical evidences. With the new implementation, the convexification algorithm is computationally efficient and is feasible to be applied in real time to detect and image antipersonnel land mines in the field.

ORBIT: Optimization by Radial Basis Function Interpolation in Trust-Regions

Stefan M. Wild, Rommel G. Regis, and Christine A. Shoemaker

SIAM J. Sci. Comput. 30, pp. 3197-3219 (23 pages) | Cited 3 times

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
We present a new derivative-free algorithm, ORBIT, for unconstrained local optimization of computationally expensive functions. A trust-region framework using interpolating Radial Basis Function (RBF) models is employed. The RBF models considered often allow ORBIT to interpolate nonlinear functions using fewer function evaluations than the polynomial models considered by present techniques. Approximation guarantees are obtained by ensuring that a subset of the interpolation points is sufficiently poised for linear interpolation. The RBF property of conditional positive definiteness yields a natural method for adding additional points. We present numerical results on test problems to motivate the use of ORBIT when only a relatively small number of expensive function evaluations are available. Results on two very different application problems, calibration of a watershed model and optimization of a PDE-based bioremediation plan, are also encouraging and support ORBIT's effectiveness on blackbox functions for which no special mathematical structure is known or available.

Multigrid One-Shot Method for State Constrained Aerodynamic Shape Optimization

Subhendu Bikash Hazra

SIAM J. Sci. Comput. 30, pp. 3220-3248 (29 pages)

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
This paper presents applications of an optimization-based multigrid method to state constrained aerodynamic shape optimization problems. The one-shot simultaneous pseudo-time-stepping method has been applied successfully to such problems, reducing the cost of computation up to 70% from that of a traditional gradient method. However, the number of optimization iterations is comparatively large in this method since we update the design parameters in each time-step of state and costate runs. In this paper we use an optimization-based multigrid strategy to reduce the total number of optimization iterations while keeping the idea of the one-shot pseudo-time-stepping method. Design examples of drag reduction with constant lift, together with geometrical constraints, for an RAE2822 airfoil are included. The total effort required for convergence of the optimization problem is between four and six times that of forward simulation runs depending on computational grid size.

Incremental Identification of Transport Coefficients in Convection-Diffusion Systems

Maka Karalashvili, Sven Groß, Adel Mhamdi, Arnold Reusken, and Wolfgang Marquardt

SIAM J. Sci. Comput. 30, pp. 3249-3269 (21 pages) | Cited 1 time

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
In this paper, an incremental approach for the identification of a model for transport coefficients in convection-diffusion systems on the basis of high-resolution measurement data is presented. The transport is represented by a convection term with known convective velocity and by a diffusion term with an unknown, generally state-dependent transport coefficient. The identification of the transport model for this transport coefficient constitutes an ill-posed nonlinear inverse problem. We present a novel decomposition approach in which this inverse problem is split into a sequence of inverse subproblems. In the first identification step of this incremental approach a source is estimated by solving an affine-linear inverse problem by means of the conjugate gradient method. In the second identification step a nonlinear inverse problem has to be solved to reconstruct a transport coefficient. A Newton-type method using the conjugate gradient method in its inner iteration is used to solve this nonlinear inverse problem of coefficient estimation. Finally, in the third identification step a transport model structure is proposed and identified on the basis of the model-free transport coefficient reconstructed in the two previous steps. The ill-posedness of each inverse problem is examined by using artificially perturbed transient simulation data and appropriate regularization techniques. The identification methodology is illustrated for a three-dimensional convection-diffusion equation which has its origin in the modeling and simulation of energy transport in a laminar wavy film flow.

Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space

T. Bui-Thanh, K. Willcox, and O. Ghattas

SIAM J. Sci. Comput. 30, pp. 3270-3288 (19 pages) | Cited 4 times

Online Publication Date: October 13, 2008

Full Text: | Download PDF

Show Abstract
A model-constrained adaptive sampling methodology is proposed for the reduction of large-scale systems with high-dimensional parametric input spaces. Our model reduction method uses a reduced basis approach, which requires the computation of high-fidelity solutions at a number of sample points throughout the parametric input space. A key challenge that must be addressed in the optimization, control, and probabilistic settings is the need for the reduced models to capture variation over this parametric input space, which, for many applications, will be of high dimension. We pose the task of determining appropriate sample points as a PDE-constrained optimization problem, which is implemented using an efficient adaptive algorithm that scales well to systems with a large number of parameters. The methodology is demonstrated using examples with parametric input spaces of dimension 11 and 21, which describe thermal analysis and design of a heat conduction fin, and compared with statistically based sampling methods. For these examples, the model-constrained adaptive sampling leads to reduced models that, for a given basis size, have error several orders of magnitude smaller than that obtained using the other methods.
Close

close