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. / Doctor of Philosophy / The study of modern science and engineering begins with descriptions of a system of mathematical equations (a model). Different models require different techniques to both accurately and effectively solve them on a computer. In this dissertation, we focus on developing novel mathematical solvers for models expressed as a system of equations, where only the initial state and the rate of change of state as a function are known. The solvers we develop can be used to both forecast the behavior of the system and to optimize its characteristics to achieve specific goals. We also build methodologies to estimate and control errors introduced by mathematical solvers in obtaining a solution for models involving multiple interacting physical, chemical, or biological phenomena. Our solvers build on state of the art in the research community by introducing new approximations that exploit the underlying mathematical structure of a model. Where it is necessary, we provide concrete mathematical proofs to validate theoretically the correctness of the approximations we introduce and correlate with follow-up experiments. We also present detailed descriptions of the procedure for implementing each mathematical solver that we develop throughout the dissertation while emphasizing on means to obtain maximal performance from the solver. We demonstrate significant performance improvements on a range of models that serve as running examples, describing chemical reactions among distinct species as they diffuse over a surface medium. Also provided are results and procedures that a curious researcher can use to advance the ideas presented in the dissertation to other types of solvers that we have not considered. Research on mathematical solvers for different mathematical models is rich and rewarding with numerous open-ended questions and is a critical component in the progress of modern science and engineering.

Accurate numerical solution of time-dependent, high-dimensional partial differential equations (PDEs) usually requires efficient numerical techniques and massive-scale parallel computing. In this thesis, we implement and evaluate discretization schemes suited for PDEs of higher dimensionality, focusing on high order of accuracy and low computational cost. Spatial discretization is particularly challenging in higher dimensions. The memory requirements for uniform grids quickly grow out of reach even on large-scale parallel computers. We utilize high-order discretization schemes and implement adaptive mesh refinement on structured hyperrectangular domains in order to reduce the required number of grid points and computational work. We allow for anisotropic (non-uniform) refinement by recursive bisection and show how to construct, manage and load balance such grids efficiently. In our numerical examples, we use finite difference schemes to discretize the PDEs. In the adaptive case we show how a stable discretization can be constructed using SBP-SAT operators. However, our adaptive mesh framework is general and other methods of discretization are viable. For integration in time, we implement exponential integrators based on the Lanczos/Arnoldi iterative schemes for eigenvalue approximations. Using adaptive time stepping and a truncated Magnus expansion, we attain high levels of accuracy in the solution at low computational cost. We further investigate alternative implementations of the Lanczos algorithm with reduced communication costs. As an example application problem, we have considered the time-dependent Schrödinger equation (TDSE). We present solvers and results for the solution of the TDSE on equidistant as well as adaptively refined Cartesian grids. / eSSENCE

The theory of partial differential equations constitutes today one of the most important topics of scientific understanding. A standard approach for solving a time-dependent partial differential equation consists in discretizing the spatial variables by finite differences or finite elements. This results in a huge system of (stiff) ordinary differential equations that has to be integrated in time. Exponential integrators constitute an interesting class of numerical methods for the time integration of stiff systems of differential equations. Their efficient implementation heavily relies on the fast computation of the action of certain matrix functions; among those, the matrix exponential is the most prominent one. In this manuscript, we go through the steps that led to the development of backward error accurate routines for computing the action of the matrix exponential.

The solution of stiff systems of Ordinary Differential Equations (ODEs), that typically arise after spatial discretization of many important evolutionary Partial Differential Equations (PDEs), constitutes a topic of wide interest in numerical analysis. A prominent way to numerically integrate such systems involves using exponential integrators. In general, these kinds of schemes do not require the solution of (non)linear systems but rather the action of the matrix exponential and of some specific exponential-like functions (known in the literature as φ-functions). In this PhD thesis we aim at presenting efficient tensor-based tools to approximate such actions, both from a theoretical and from a practical point of view, when the problem has an underlying Kronecker sum structure. Moreover, we investigate the application of exponential integrators to compute numerical solutions of important equations in various fields, such as plasma physics, mean-field optimal control and computational chemistry. In any case, we provide several numerical examples and we perform extensive simulations, eventually exploiting modern hardware architectures such as multi-core Central Processing Units (CPUs) and Graphic Processing Units (GPUs). The results globally show the effectiveness and the superiority of the different approaches proposed.

The time-dependent Schrödinger equation (TDSE) models the quantum nature of molecular processes.  Numerical simulations based on the TDSE help in understanding and predicting the outcome of chemical reactions. This thesis is dedicated to the derivation and analysis of efficient and reliable simulation tools for the TDSE, with a particular focus on models for the interaction of molecules with time-dependent electromagnetic fields. Various time propagators are compared for this setting and an efficient fourth-order commutator-free Magnus-Lanczos propagator is derived. For the Lanczos method, several communication-reducing variants are studied for an implementation on clusters of multi-core processors. Global error estimation for the Magnus propagator is devised using a posteriori error estimation theory. In doing so, the self-adjointness of the linear Schrödinger equation is exploited to avoid solving an adjoint equation. Efficiency and effectiveness of the estimate are demonstrated for both bounded and unbounded states. The temporal approximation is combined with adaptive spectral elements in space. Lagrange elements based on Gauss-Lobatto nodes are employed to avoid nondiagonal mass matrices and ill-conditioning at high order. A matrix-free implementation for the evaluation of the spectral element operators is presented. The framework uses hybrid parallelism and enables significant computational speed-up as well as the solution of larger problems compared to traditional implementations relying on sparse matrices. As an alternative to grid-based methods, radial basis functions in a Galerkin setting are proposed and analyzed. It is found that considerably higher accuracy can be obtained with the same number of basis functions compared to the Fourier method. Another direction of research presented in this thesis is a new algorithm for quantum optimal control: The field is optimized in the frequency domain where the dimensionality of the optimization problem can drastically be reduced. In this way, it becomes feasible to use a quasi-Newton method to solve the problem. / eSSENCE