Computational Science Laboratory

Permanent URI for this community

https://hdl.handle.net/10919/102074

The mission of the Computational Science Laboratory (CSL) is to develop innovative computational solutions for complex real-world problems, and to foster a productive research and education environment emphasizing collaboration and innovation.

Browse

Now showing 1 - 20 of 131

Accelerating Atmospheric Modeling Through Emerging Multi-core Technologies
Linford, John Christian (Virginia Tech, 2010-05-05)
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.
Achieving Very High Order for Implicit Explicit Time Stepping: Extrapolation Methods
Constantinescu, Emil M.; Sandu, Adrian (Department of Computer Science, Virginia Polytechnic Institute & State University, 2008-07-01)
In this paper we construct extrapolated implicit-explicit time stepping methods that allow to efficiently solve problems with both stiff and non-stiff components. The proposed methods can provide very high order discretizations of ODEs, index-1 DAEs, and PDEs in the method of lines framework. These methods are simple to construct, easy to implement and parallelize. We establish the existence of perturbed asymptotic expansions of global errors, explain the convergence orders of these methods, and explore their linear stability properties. Numerical results with stiff ODEs, DAEs, and PDEs illustrate the theoretical findings and the potential of these methods to solve multiphysics multiscale problems.
Adaptive Numerical Methods for Large Scale Simulations and Data Assimilation
Constantinescu, Emil Mihai (Virginia Tech, 2008-05-26)
Numerical simulation is necessary to understand natural phenomena, make assessments and predictions in various research and engineering fields, develop new technologies, etc. New algorithms are needed to take advantage of the increasing computational resources and utilize the emerging hardware and software infrastructure with maximum efficiency. Adaptive numerical discretization methods can accommodate problems with various physical, scale, and dynamic features by adjusting the resolution, order, and the type of method used to solve them. In applications that simulate real systems, the numerical accuracy of the solution is typically just one of the challenges. Measurements can be included in the simulation to constrain the numerical solution through a process called data assimilation in order to anchor the simulation in reality. In this thesis we investigate adaptive discretization methods and data assimilation approaches for large-scale numerical simulations. We develop and investigate novel multirate and implicit-explicit methods that are appropriate for multiscale and multiphysics numerical discretizations. We construct and explore data assimilation approaches for, but not restricted to, atmospheric chemistry applications. A generic approach for describing the structure of the uncertainty in initial conditions that can be applied to the most popular data assimilation approaches is also presented. We show that adaptive numerical methods can effectively address the discretization of large-scale problems. Data assimilation complements the adaptive numerical methods by correcting the numerical solution with real measurements. Test problems and large-scale numerical experiments validate the theoretical findings. Synergistic approaches that use adaptive numerical methods within a data assimilation framework need to be investigated in the future.
Adjoint based solution and uncertainty quantification techniques for variational inverse problems
Hebbur Venkata Subba Rao, Vishwas (Virginia Tech, 2015-09-25)
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 second order adjoints and provides an a-posteriori error estimate for a quantity of interest defined on the inverse solution (i.e., an aspect of the inverse solution). We illustrate the estimation algorithm on a shallow water model and on the Weather Research and Forecast model. Presence of outliers in measurement data is common, and this negatively impacts the solution to variational inverse problems. The traditional approach, where the inverse problem is formulated as a minimization problem in $L_2$ norm, is especially sensitive to large data errors. To alleviate the impact of data outliers we propose to use robust norms such as the $L_1$ and Huber norm in data assimilation. This work develops a systematic mathematical framework to perform three and four dimensional variational data assimilation using $L_1$ and Huber norms. The power of this approach is demonstrated by solving data assimilation problems where measurements contain outliers.
Adjoint-based space-time adaptive solution algorithms for sensitivity analysis and inverse problems
Alexe, Mihai (Virginia Tech, 2011-03-18)
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.
Adjoint-Matching Neural Network Surrogates for Fast 4D-Var Data Assimilation
Chennault, Austin; Popov, Andrey A.; Subrahmanya, Amit N.; Cooper, Rachel; Karpatne, Anuj; Sandu, Adrian (2021-11-16)
The data assimilation procedures used in many operational numerical weather forecasting systems are based around variants of the 4D-Var algorithm. The cost of solving the 4D-Var problem is dominated by the cost of forward and adjoint evaluations of the physical model. This motivates their substitution by fast, approximate surrogate models. Neural networks offer a promising approach for the data-driven creation of surrogate models. The accuracy of the surrogate 4D-Var problem’s solution has been shown to depend explicitly on accurate modeling of the forward and adjoint for other surrogate modeling approaches and in the general nonlinear setting. We formulate and analyze several approaches to incorporating derivative information into the construction of neural network surrogates. The resulting networks are tested on out of training set data and in a sequential data assimilation setting on the Lorenz-63 system. Two methods demonstrate superior performance when compared with a surrogate network trained without adjoint information, showing the benefit of incorporating adjoint information into the training process.
Advanced Sampling Methods for Solving Large-Scale Inverse Problems
Attia, Ahmed Mohamed Mohamed (Virginia Tech, 2016-09-19)
Ensemble and variational techniques have gained wide popularity as the two main approaches for solving data assimilation and inverse problems. The majority of the methods in these two approaches are derived (at least implicitly) under the assumption that the underlying probability distributions are Gaussian. It is well accepted, however, that the Gaussianity assumption is too restrictive when applied to large nonlinear models, nonlinear observation operators, and large levels of uncertainty. This work develops a family of fully non-Gaussian data assimilation algorithms that work by directly sampling the posterior distribution. The sampling strategy is based on a Hybrid/Hamiltonian Monte Carlo (HMC) approach that can handle non-normal probability distributions. The first algorithm proposed in this work is the "HMC sampling filter", an ensemble-based data assimilation algorithm for solving the sequential filtering problem. Unlike traditional ensemble-based filters, such as the ensemble Kalman filter and the maximum likelihood ensemble filter, the proposed sampling filter naturally accommodates non-Gaussian errors and nonlinear model dynamics, as well as nonlinear observations. To test the capabilities of the HMC sampling filter numerical experiments are carried out using the Lorenz-96 model and observation operators with different levels of nonlinearity and differentiability. The filter is also tested with shallow water model on the sphere with linear observation operator. Numerical results show that the sampling filter performs well even in highly nonlinear situations where the traditional filters diverge. Next, the HMC sampling approach is extended to the four-dimensional case, where several observations are assimilated simultaneously, resulting in the second member of the proposed family of algorithms. The new algorithm, named "HMC sampling smoother", is an ensemble-based smoother for four-dimensional data assimilation that works by sampling from the posterior probability density of the solution at the initial time. The sampling smoother naturally accommodates non-Gaussian errors and nonlinear model dynamics and observation operators, and provides a full description of the posterior distribution. Numerical experiments for this algorithm are carried out using a shallow water model on the sphere with observation operators of different levels of nonlinearity. The numerical results demonstrate the advantages of the proposed method compared to the traditional variational and ensemble-based smoothing methods. The HMC sampling smoother, in its original formulation, is computationally expensive due to the innate requirement of running the forward and adjoint models repeatedly. The proposed family of algorithms proceeds by developing computationally efficient versions of the HMC sampling smoother based on reduced-order approximations of the underlying model dynamics. The reduced-order HMC sampling smoothers, developed as extensions to the original HMC smoother, are tested numerically using the shallow-water equations model in Cartesian coordinates. The results reveal that the reduced-order versions of the smoother are capable of accurately capturing the posterior probability density, while being significantly faster than the original full order formulation. In the presence of nonlinear model dynamics, nonlinear observation operator, or non-Gaussian errors, the prior distribution in the sequential data assimilation framework is not analytically tractable. In the original formulation of the HMC sampling filter, the prior distribution is approximated by a Gaussian distribution whose parameters are inferred from the ensemble of forecasts. The Gaussian prior assumption in the original HMC filter is relaxed. Specifically, a clustering step is introduced after the forecast phase of the filter, and the prior density function is estimated by fitting a Gaussian Mixture Model (GMM) to the prior ensemble. The base filter developed following this strategy is named cluster HMC sampling filter (ClHMC ). A multi-chain version of the ClHMC filter, namely MC-ClHMC , is also proposed to guarantee that samples are taken from the vicinities of all probability modes of the formulated posterior. These methodologies are tested using a quasi-geostrophic (QG) model with double-gyre wind forcing and bi-harmonic friction. Numerical results demonstrate the usefulness of using GMMs to relax the Gaussian prior assumption in the HMC filtering paradigm. To provide a unified platform for data assimilation research, a flexible and a highly-extensible testing suite, named DATeS , is developed and described in this work. The core of DATeS is implemented in Python to enable for Object-Oriented capabilities. The main components, such as the models, the data assimilation algorithms, the linear algebra solvers, and the time discretization routines are independent of each other, such as to offer maximum flexibility to configure data assimilation studies.
Advanced Time Integration Methods with Applications to Simulation, Inverse Problems, and Uncertainty Quantification
Narayanamurthi, Mahesh (Virginia Tech, 2020-01-29)
Simulation and optimization of complex physical systems are an integral part of modern science and engineering. The systems of interest in many fields have a multiphysics nature, with complex interactions between physical, chemical and in some cases even biological processes. This dissertation seeks to advance forward and adjoint numerical time integration methodologies for the simulation and optimization of semi-discretized multiphysics partial differential equations (PDEs), and to estimate and control numerical errors via a goal-oriented a posteriori error framework. We extend exponential propagation iterative methods of Runge-Kutta type (EPIRK) by [Tokman, JCP 2011], to build EPIRK-W and EPIRK-K time integration methods that admit approximate Jacobians in the matrix-exponential like operations. EPIRK-W methods extend the W-method theory by [Steihaug and Wofbrandt, Math. Comp. 1979] to preserve their order of accuracy under arbitrary Jacobian approximations. EPIRK-K methods extend the theory of K-methods by [Tranquilli and Sandu, JCP 2014] to EPIRK and use a Krylov-subspace based approximation of Jacobians to gain computational efficiency. New families of partitioned exponential methods for multiphysics problems are developed using the classical order condition theory via particular variants of T-trees and corresponding B-series. The new partitioned methods are found to perform better than traditional unpartitioned exponential methods for some problems in mild-medium stiffness regimes. Subsequently, partitioned stiff exponential Runge-Kutta (PEXPRK) methods -- that extend stiffly accurate exponential Runge-Kutta methods from [Hochbruck and Ostermann, SINUM 2005] to a multiphysics context -- are constructed and analyzed. PEXPRK methods show full convergence under various splittings of a diffusion-reaction system. We address the problem of estimation of numerical errors in a multiphysics discretization by developing a goal-oriented a posteriori error framework. Discrete adjoints of GARK methods are derived from their forward formulation [Sandu and Guenther, SINUM 2015]. Based on these, we build a posteriori estimators for both spatial and temporal discretization errors. We validate the estimators on a number of reaction-diffusion systems and use it to simultaneously refine spatial and temporal grids.
Alternating directions implicit integration in a general linear method framework
Sarshar, Arash; Roberts, Steven; Sandu, Adrian (Elsevier, 2021-05-15)
Alternating Directions Implicit (ADI) integration is an operator splitting approach to solve parabolic and elliptic partial differential equations in multiple dimensions based on solving sequentially a set of related one-dimensional equations. Classical ADI methods have order at most two, due to the splitting errors. Moreover, when the time discretization of stiff one-dimensional problems is based on Runge–Kutta schemes, additional order reduction may occur. This work proposes a new ADI approach based on the partitioned General Linear Methods framework. This approach allows the construction of high order ADI methods. Due to their high stage order, the proposed methods can alleviate the order reduction phenomenon seen with other schemes. Numerical experiments are shown to provide further insight into the accuracy, stability, and applicability of these new methods.
Application of approximate matrix factorization to high order linearly implicit Runge-Kutta methods
Zhang, H.; Sandu, Adrian; Tranquilli, Paul (Elsevier, 2015-10-01)
Approximate Exponential Algorithms to Solve the Chemical Master Equation
Mooasvi, A.; Sandu, Adrian (Vilnius Gediminas Tech Univ, 2015-05-04)
Augmented Neural Network Surrogate Models for Polynomial Chaos Expansions and Reduced Order Modeling
Cooper, Rachel Gray (Virginia Tech, 2021-05-20)
Mathematical models describing real world processes are becoming increasingly complex to better match the dynamics of the true system. While this is a positive step towards more complete knowledge of our world, numerical evaluations of these models become increasingly computationally inefficient, requiring increased resources or time to evaluate. This has led to the need for simplified surrogates to these complex mathematical models. A growing surrogate modeling solution is with the usage of neural networks. Neural networks (NN) are known to generalize an approximation across a diverse dataset and minimize the solution along complex nonlinear boundaries. Additionally, these surrogate models can be found using only incomplete knowledge of the true dynamics. However, NN surrogates often suffer from a lack of interpretability, where the decisions made in the training process are not fully understood, and the roles of individual neurons are not well defined. We present two solutions towards this lack of interpretability. The first focuses on mimicking polynomial chaos (PC) modeling techniques, modifying the structure of a NN to produce polynomial approximations of the underlying dynamics. This methodology allows for an extractable meaning from the network and results in improvement in accuracy over traditional PC methods. Secondly, we examine the construction of a reduced order modeling scheme using NN autoencoders, guiding the decisions of the training process to better match the real dynamics. This guiding process is performed via a physics-informed (PI) penalty, resulting in a speed-up in training convergence, but still results in poor performance compared to traditional schemes.
Autoregressive Models of Background Errors for Chemical Data Assimilation
Constantinescu, Emil M.; Chai, Tianfeng; Sandu, Adrian; Carmichael, Gregory R. (Department of Computer Science, Virginia Polytechnic Institute & State University, 2006-10-01)
The task of providing an optimal analysis of the state of the atmosphere requires the development of dynamic data-driven systems that efficiently integrate the observational data and the models. Data assimilation (DA) is the process of adjusting the states or parameters of a model in such a way that its outcome (prediction) is close, in some distance metric, to observed (real) states. It is widely accepted that a key ingredient of successful data assimilation is a realistic estimation of the background error distribution. This paper introduces a new method for estimating the background errors which are modeled using autoregressive processes. The proposed approach is computationally inexpensive and captures the error correlations along the flow lines.
A Bayesian approach to multivariate adaptive localization in ensemble-based data assimilation with time-dependent extensions
Popov, Andrey A.; Sandu, Adrian (Copernicus Publications, 2019-06-14)
Ever since its inception, the ensemble Kalman filter (EnKF) has elicited many heuristic approaches that sought to improve it. One such method is covariance localization, which alleviates spurious correlations due to finite ensemble sizes by using relevant spatial correlation information. Adaptive localization techniques account for how correlations change in time and space, in order to obtain improved covariance estimates. This work develops a Bayesian approach to adaptive Schur-product localization for the deterministic ensemble Kalman filter (DEnKF) and extends it to support multiple radii of influence. We test the proposed adaptive localization using the toy Lorenz’96 problem and a more realistic 1.5-layer quasi-geostrophic model. Results with the toy problem show that the multivariate approach informs us that strongly observed variables can tolerate larger localization radii. The univariate approach leads to markedly improved filter performance for the realistic geophysical model, with a reduction in error by as much as 33 %.
Chemical Data Assimilation-An Overview
Sandu, Adrian; Chai, Tianfeng (MDPI, 2011-08-29)
Chemical data assimilation is the process by which models use measurements to produce an optimal representation of the chemical composition of the atmosphere. Leveraging advances in algorithms and increases in the available computational power, the integration of numerical predictions and observations has started to play an important role in air quality modeling. This paper gives an overview of several methodologies used in chemical data assimilation. We discuss the Bayesian framework for developing data assimilation systems, the suboptimal and the ensemble Kalman filter approaches, the optimal interpolation (OI), and the three and four dimensional variational methods. Examples of assimilation real observations with CMAQ model are presented.
Chemical Mechanism Solvers in Air Quality Models
Zhang, Hong; Linford, John C.; Sandu, Adrian; Sander, Rolf (MDPI, 2011-09-01)
The solution of chemical kinetics is one of the most computationally intensive tasks in atmospheric chemical transport simulations. Due to the stiff nature of the system, implicit time stepping algorithms which repeatedly solve linear systems of equations are necessary. This paper reviews the issues and challenges associated with the construction of efficient chemical solvers, discusses several families of algorithms, presents strategies for increasing computational efficiency, and gives insight into implementing chemical solvers on accelerated computer architectures.
A class of generalized additive Runge-Kutta methods
Sandu, Adrian; Guenther, Michael (2013-10-22)
This work generalizes the additively partitioned Runge-Kutta methods by allowing for different stage values as arguments of different components of the right hand side. An order conditions theory is developed for the new family of generalized additive methods, and stability and monotonicity investigations are carried out. The paper discusses the construction and properties of implicit-explicit and implicit-implicit,methods in the new framework. The new family, named GARK, introduces additional flexibility when compared to traditional partitioned Runge-Kutta methods, and therefore offers additional opportunities for the development of flexible solvers for systems with multiple scales, or driven by multiple physical processes.
A class of implicit-explicit two-step Runge-Kutta methods
Zharovski, Evgeniy; Sandu, Adrian (Department of Computer Science, Virginia Polytechnic Institute & State University, 2012-02-01)
This work develops implicit-explicit time integrators based on two-step Runge-Kutta methods. The class of schemes of interest is characterized by linear invariant preservation and high stage orders. Theoretical consistency and stability analyses are performed to reveal the properties of these methods. The new framework offers extreme flexibility in the construction of partitioned integrators, since no coupling conditions are necessary. Moreover, the methods are not plagued by severe order reduction, due to their high stage orders. Two practical schemes of orders four and six are constructed, and are used to solve several test problems. Numerical results confirm the theoretical findings.
Cluster Sampling Filters for Non-Gaussian Data Assimilation
Attia, A.; Moosavi, Azam; Sandu, Adrian (2016-08-19)
This paper presents a fully non-Gaussian version of the Hamiltonian Monte Carlo (HMC) sampling filter. The Gaussian prior assumption in the original HMC filter is relaxed. Specifically, a clustering step is introduced after the forecast phase of the filter, and the prior density function is estimated by fitting a Gaussian Mixture Model (GMM) to the prior ensemble. Using the data likelihood function, the posterior density is then formulated as a mixture density, and is sampled using a HMC approach (or any other scheme capable of sampling multimodal densities in high-dimensional subspaces). The main filter developed herein is named "cluster HMC sampling filter" (ClHMC). A multi-chain version of the ClHMC filter, namely MC-ClHMC is also proposed to guarantee that samples are taken from the vicinities of all probability modes of the formulated posterior. The new methodologies are tested using a quasi-geostrophic (QG) model with double-gyre wind forcing and bi-harmonic friction. Numerical results demonstrate the usefulness of using GMMs to relax the Gaussian prior assumption in the HMC filtering paradigm.
Cluster Sampling Filters for Non-Gaussian Data Assimilation
Attia, Ahmed; Moosavi, Azam; Sandu, Adrian (MDPI, 2018-05-31)
This paper presents a fully non-Gaussian filter for sequential data assimilation. The filter is named the “cluster sampling filter”, and works by directly sampling the posterior distribution following a Markov Chain Monte-Carlo (MCMC) approach, while the prior distribution is approximated using a Gaussian Mixture Model (GMM). Specifically, a clustering step is introduced after the forecast phase of the filter, and the prior density function is estimated by fitting a GMM to the prior ensemble. Using the data likelihood function, the posterior density is then formulated as a mixture density, and is sampled following an MCMC approach. Four versions of the proposed filter, namely C ℓ MCMC , C ℓ HMC , MC- C ℓ HMC , and MC- C ℓ HMC are presented. C ℓ MCMC uses a Gaussian proposal density to sample the posterior, and C ℓ HMC is an extension to the Hamiltonian Monte-Carlo (HMC) sampling filter. MC- C ℓ MCMC and MC- C ℓ HMC are multi-chain versions of the cluster sampling filters C ℓ MCMC and C ℓ HMC respectively. The multi-chain versions are proposed to guarantee that samples are taken from the vicinities of all probability modes of the formulated posterior. The new methodologies are tested using a simple one-dimensional example, and a quasi-geostrophic (QG) model with double-gyre wind forcing and bi-harmonic friction. Numerical results demonstrate the usefulness of using GMMs to relax the Gaussian prior assumption especially in the HMC filtering paradigm.

Browse

Browsing Computational Science Laboratory by Title

Results Per Page

Sort Options