Search

Mathematical models involving ordinary differential equations (ODEs) play a critical role in scientific and engineering applications. Advances in computing hardware and numerical methods have allowed these models to become larger and more sophisticated. Increasingly, problems can be described as multiphysics and multiscale as they combine several different physical processes with different characteristics. If just one part of an ODE is stiff, nonlinear, chaotic, or rapidly-evolving, this can force an expensive method or a small timestep to be used. A method which applies a discretization and timestep uniformly across a multiphysics problem poorly utilizes computational resources and can be prohibitively expensive. The focus of this dissertation is on "multimethods" which apply different methods to different partitions of an ODE. Well-designed multimethods can drastically reduce the computation costs by matching methods to the individual characteristics of each partition while making minimal concessions to stability and accuracy. However, they are not without their limitations. High order methods are difficult to derive and may suffer from order reduction. Also, the stability of multimethods is difficult to characterize and analyze. The goals of this work are to develop new, practical multimethods and to address these issues. First, new implicit multirate Runge–Kutta methods are analyzed with a special focus on stability. This is extended into implicit multirate infinitesimal methods. We introduce approaches for constructing implicit-explicit methods based on Runge–Kutta and general linear methods. Finally, some unique applications of multimethods are considered including using surrogate models to accelerate Runge–Kutta methods and eliminating order reduction on linear ODEs with time-dependent forcing.

Mathematical modeling of physical processes often leads to systems of differential and algebraic equations involving quantities of interest. A computer model created based on these equations can be numerically integrated to predict future states of the system and its evolution in time. This thesis investigates current methods in numerical time-stepping schemes, identifying a number of important features needed to speed up and increase the accuracy of the solutions. The focus is on developing new methods suitable for large-scale applications with multiple physical processes, potentially with significant differences in their time-scales. Various families of new methods are introduced with special attention to multirating, low computational cost implicitness, high order of convergence, and robustness. For each family, the order condition theory is discussed and a number of examples are derived. The accuracy and stability of the methods are investigated using standard analysis techniques and numerical experiments are performed to verify the abilities of the new methods.

The solution of initial value problems is a fundamental component of many scientific simulations of physical phenomena. In many cases these initial value problems arise from a method of lines approach to solving partial differential equations, resulting in very large systems of equations that require the use of numerical time integration methods to solve. Many problems of scientific interest exhibit stiff behavior for which implicit methods are favorable, however standard implicit methods are computationally expensive. They require the solution of one or more large nonlinear systems at each timestep, which can be impractical to solve exactly and can behave poorly when solved approximately. The recently introduced ``lightly-implicit'' K-methods seek to avoid this issue by directly coupling the time integration methods with a Krylov based approximation of linear system solutions, treating a portion of the problem implicitly and the remainder explicitly. This work seeks to further two primary objectives: evaluation of these K-methods in large-scale parallel applications, and development of new linearly implicit methods for contexts where improvements can be made. To this end, Rosenbrock-Krylov methods, the first K-methods, are examined in a scalability study, and two new families of time integration methods are introduced: biorthogonal Rosenbrock-Krylov methods, and linearly implicit multistep methods. For the scalability evaluation of Rosenbrock-Krylov methods, two parallel contexts are considered: a GPU accelerated model and a distributed MPI parallel model. In both cases, the most significant performance bottleneck is the need for many vector dot products, which require costly parallel reduce operations. Biorthogonal Rosenbrock-Krylov methods are an extension of the original Rosenbrock-Krylov methods which replace the Arnoldi iteration used to produce the Krylov approximation with Lanczos biorthogonalization, which requires fewer vector dot products, leading to lowe

Simulations and modeling of large-scale systems are vital to understanding real world phenomena. However, even advanced numerical models can only approximate the true physics. The discrepancy between model results and nature can be attributed to different sources of uncertainty including the parameters of the model, input data, or some missing physics that is not included in the model due to a lack of knowledge or high computational costs. Uncertainty reduction approaches seek to improve the model accuracy by decreasing the overall uncertainties in models. Aiming to contribute to this area, this study explores uncertainty quantification and reduction approaches for complex physical problems. This study proposes several novel probabilistic and statistical approaches for identifying the sources of uncertainty, modeling the errors, and reducing uncertainty to improve the model predictions for large-scale simulations. We explore different computational models. The first class of models studied herein are inherently stochastic, and numerical approximations suffer from stability and accuracy issues. The second class of models are partial differential equations, which capture the laws of mathematical physics; however, they only approximate a more complex reality, and have uncertainties due to missing dynamics which is not captured by the models. The third class are low-fidelity models, which are fast approximations of very expensive high-fidelity models. The reduced-order models have uncertainty due to loss of information in the dimension reduction process. We also consider uncertainty analysis in the data assimilation framework, specifically for ensemble based methods where the effect of sampling errors is alleviated by localization. Finally, we study the uncertainty in numerical weather prediction models coming from approximate descriptions of physical processes.

Enhance cybersecurity learning experiences of students at Virginia Tech’s large engineering program

Ensemble-based methods have gained widespread popularity in the field of data assimilation. An ensemble of model realizations encapsulates information about the error correlations driven by the physics and the dynamics of the numerical model. This information can be used to obtain improved estimates of the state of non-linear dynamical systems such as the atmosphere and/or the ocean. This work develops efficient ensemble-based methods for data assimilation. A major bottleneck in ensemble Kalman filter (EnKF) implementations is the solution of a linear system at each analysis step. To alleviate it an EnKF implementation based on an iterative Sherman Morrison formula is proposed. The rank deficiency of the ensemble covariance matrix is exploited in order to efficiently compute the analysis increments during the assimilation process. The computational effort of the proposed method is comparable to those of the best EnKF implementations found in the current literature. The stability analysis of the new algorithm is theoretically proven based on the positiveness of the data error covariance matrix. In order to improve the background error covariance matrices in ensemble-based data assimilation we explore the use of shrinkage covariance matrix estimators from ensembles. The resulting filter has attractive features in terms of both memory usage and computational complexity. Numerical results show that it performs better that traditional EnKF formulations. In geophysical applications the correlations between errors corresponding to distant model components decreases rapidly with the distance. We propose a new and efficient implementation of the EnKF based on a modified Cholesky decomposition for inverse covariance matrix estimation. This approach exploits the conditional independence of background errors between distant model components with regard to a predefined radius of influence. Consequently, sparse estimators of the inverse background error covariance matrix can be o

Many scientific and engineering applications require the solution of large systems of initial value problems arising from method of lines discretization of partial differential equations. For systems with widely varying time scales, or with complex physical dynamics, implicit time integration schemes are preferred due to their superior stability properties. However, for very large systems accurate solution of the implicit terms can be impractical. For this reason approximations are widely used in the implementation of such methods. The primary focus of this work is on the development of novel ``lightly-implicit'' time integration methodologies. These methods consider the time integration and the solution of the implicit terms as a single computational process. We propose several classes of lightly-implicit methods that can be constructed to allow for different, specific approximations. Rosenbrock-Krylov and exponential-Krylov methods are designed to permit low accuracy Krylov based approximations of the implicit terms, while maintaining full order of convergence. These methods are matrix free, have low memory requirements, and are particularly well suited to parallel architectures. Linear stability analysis of K-methods is leveraged to construct implementation improvements for both Rosenbrock-Krylov and exponential-Krylov methods. Linearly-implicit Runge-Kutta-W methods are designed to permit arbitrary, time dependent, and stage varying approximations of the linear stiff dynamics of the initial value problem. The methods presented here are constructed with approximate matrix factorization in mind, though the framework is flexible and can be extended to many other approximations. The flexibility of lightly-implicit methods, and their ability to leverage computationally favorable approximations makes them an ideal alternative to standard explicit and implicit schemes for large parallel applications.

Variational inverse problems integrate computational simulations of physical phenomena with physical measurements in an informational feedback control system. Control parameters of the computational model are optimized such that the simulation results fit the physical measurements.The solution procedure is computationally expensive since it involves running the simulation computer model (the emph{forward model}) and the associated emph {adjoint model} multiple times. In practice, our knowledge of the underlying physics is incomplete and hence the associated computer model is laden with emph {model errors}. Similarly, it is not possible to measure the physical quantities exactly and hence the measurements are associated with emph {data errors}. The errors in data and model adversely affect the inference solutions. This work develops methods to address the challenges posed by the computational costs and by the impact of data and model errors in solving variational inverse problems. Variational inverse problems of interest here are formulated as optimization problems constrained by partial differential equations (PDEs). The solution process requires multiple evaluations of the constraints, therefore multiple solutions of the associated PDE. To alleviate the computational costs we develop a parallel in time discretization algorithm based on a nonlinear optimization approach. Like in the emph{parareal} approach, the time interval is partitioned into subintervals, and local time integrations are carried out in parallel. Solution continuity equations across interval boundaries are added as constraints. All the computational steps - forward solutions, gradients, and Hessian-vector products - involve only ideally parallel computations and therefore are highly scalable. This work develops a systematic mathematical framework to compute the impact of data and model errors on the solution to the variational inverse problems. The computational algorithm makes use of first and sec

Many fields in science and engineering require large-scale numerical simulations of complex systems described by differential equations. These systems are typically multi-physics (they are driven by multiple interacting physical processes) and multiscale (the dynamics takes place on vastly different spatial and temporal scales). Numerical solution of such systems is highly challenging due to the dimension of the resulting discrete problem, and to the complexity that comes from incorporating multiple interacting components with different characteristics. The main contributions of this dissertation are the creation of new families of time integration methods for multiscale and multiphysics simulations, and the development of industrial-strengh tools for sensitivity analysis. This work develops novel implicit-explicit (IMEX) general linear time integration methods for multiphysics and multiscale simulations typically involving both stiff and non-stiff components. In an IMEX approach, one uses an implicit scheme for the stiff components and an explicit scheme for the non-stiff components such that the combined method has the desired stability and accuracy properties. Practical schemes with favorable properties, such as maximized stability, high efficiency, and no order reduction, are constructed and applied in extensive numerical experiments to validate the theoretical findings and to demonstrate their advantages. Approximate matrix factorization (AMF) technique exploits the structure of the Jacobian of the implicit parts, which may lead to further efficiency improvement of IMEX schemes. We have explored the application of AMF within some high order IMEX Runge-Kutta schemes in order to achieve high efficiency. Sensitivity analysis gives quantitative information about the changes in a dynamical model outputs caused by caused by small changes in the model inputs. This information is crucial for data assimilation, model-constrained optimization, inverse problems, and uncer

A deep scientific understanding of complex physical systems, such as the atmosphere, can be achieved neither by direct measurements nor by numerical simulations alone. Data assimilation is a rigorous procedure to fuse information from a priori knowledge of the system state, the physical laws governing the evolution of the system, and real measurements, all with associated error statistics. Data assimilation produces best (a posteriori) estimates of model states and parameter values, and results in considerably improved computer simulations. The acquisition and use of observations in data assimilation raises several important scientific questions related to optimal sensor network design, quantification of data impact, pruning redundant data, and identifying the most beneficial additional observations. These questions originate in operational data assimilation practice, and have started to attract considerable interest in the recent past. This dissertation advances the state of knowledge in four dimensional variational (4D-Var) - data assimilation by developing, implementing, and validating a novel computational framework for estimating observation impact and for optimizing sensor networks. The framework builds on the powerful methodologies of second-order adjoint modeling and the 4D-Var sensitivity equations. Efficient computational approaches for quantifying the observation impact include matrix free linear algebra algorithms and low-rank approximations of the sensitivities to observations. The sensor network configuration problem is formulated as a meta-optimization problem. Best values for parameters such as sensor location are obtained by optimizing a performance criterion, subject to the constraint posed by the 4D-Var optimization. Tractable computational solutions to this "optimization-constrained" optimization problem are provided. The results of this work can be directly applied to the deployment of intelligent sensors and adaptive observations, as well as to

High performance computing (HPC) is an area of computer science and engineering that has always evolved rapidly---sometimes leading and sometimes riding succeeding waves of technical innovation. While HPC application developers and users have continued to benefit from the increasing power of these high-end resources, the increasing complexity of HPC execution environments will require more and more reliance on runtime systems. Parallelism, load-balancing, power, fault-tolerance, and hardware heterogeneity are just a few of the emerging dominant issues that require runtime solutions. In this talk I will briefly describe some of the motivations and trends in runtime systems for HPC. I will then describe two recent projects we have worked on at Virginia Tech. The first, ReSHAPE, is a runtime system that allows the number of nodes assigned to job running on a cluster to be changed at run time. Experimental results from a prototype implementation of ReSHAPE illustrate the potential of "malleable" jobs for improving overall cluster utilization and reducing turn-around time for individual jobs. The second project, Samhita, is a distributed shared memory (DSM) execution environment, which allow programs based on the widely used Pthreads library for shared memory thread parallelism to be easily ported to a distributed memory (cluster) platform. Samhita not only allows a wide range of parallel codes to be ported to a new context, but its design reduces the problem of DSM to a cache management problem, with corresponding opportunities for exploiting locality at runtime. Bio: Cal Ribbens is Associate Professor and Associate Department Head for Undergraduate Studies in the Department of Computer Science at Virginia Tech. He received a B.S. in Mathematics from Calvin College (1981) and a Ph.D. in Computer Sciences from Purdue University (1986). His research interests include parallel computation, numerical algorithms, mathematical software, and tools and environments for high per

Adaptivity in both space and time has become the norm for solving problems modeled by partial differential equations. The size of the discretized problem makes uniformly refined grids computationally prohibitive. Adaptive refinement of meshes and time steps allows to capture the phenomena of interest while keeping the cost of a simulation tractable on the current hardware. Many fields in science and engineering require the solution of inverse problems where parameters for a given model are estimated based on available measurement information. In contrast to forward (regular) simulations, inverse problems have not extensively benefited from the adaptive solver technology. Previous research in inverse problems has focused mainly on the continuous approach to calculate sensitivities, and has typically employed fixed time and space meshes in the solution process. Inverse problem solvers that make exclusive use of uniform or static meshes avoid complications such as the differentiation of mesh motion equations, or inconsistencies in the sensitivity equations between subdomains with different refinement levels. However, this comes at the cost of low computational efficiency. More efficient computations are possible through judicious use of adaptive mesh refinement, adaptive time steps, and the discrete adjoint method. This dissertation develops a complete framework for fully discrete adjoint sensitivity analysis and inverse problem solutions, in the context of time dependent, adaptive mesh, and adaptive step models. The discrete framework addresses all the necessary ingredients of a state–of–the–art adaptive inverse solution algorithm: adaptive mesh and time step refinement, solution grid transfer operators, a priori and a posteriori error analysis and estimation, and discrete adjoints for sensitivity analysis of flux–limited numerical algorithms.

The new generations of multi-core chipset architectures achieve unprecedented levels of computational power while respecting physical and economical constraints. The cost of this power is bewildering program complexity. Atmospheric modeling is a grand-challenge problem that could make good use of these architectures if they were more accessible to the average programmer. To that end, software tools and programming methodologies that greatly simplify the acceleration of atmospheric modeling and simulation with emerging multi-core technologies are developed. A general model is developed to simulate atmospheric chemical transport and atmospheric chemical kinetics. The Cell Broadband Engine Architecture (CBEA), General Purpose Graphics Processing Units (GPGPUs), and homogeneous multi-core processors (e.g. Intel Quad-core Xeon) are introduced. These architectures are used in case studies of transport modeling and kinetics modeling and demonstrate per-kernel speedups as high as 40x. A general analysis and code generation tool for chemical kinetics called "KPPA" is developed. KPPA generates highly tuned C, Fortran, or Matlab code that uses every layer of heterogeneous parallelism in the CBEA, GPGPU, and homogeneous multi-core architectures. A scalable method for simulating chemical transport is also developed. The Weather Research and Forecasting Model with Chemistry (WRF-Chem) is accelerated with these methods with good results: real forecasts of air quality are generated for the Eastern United States 65% faster than the state-of-the-art models.