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 lower overall cost for stiff problems. Linearly implicit multistep methods are a new family of implicit multistep methods that require only a single linear solve per timestep; the family includes W- and K-method variants, which admit arbitrary or Krylov based approximations of the problem Jacobian while maintaining the order of accuracy. This property allows for a wide range of implementation optimizations. Finally, all the new methods proposed herein are implemented efficiently in the MATLODE package, a Matlab ODE solver and sensitivity analysis toolbox, to make them available to the community at large. / Doctor of Philosophy / Differential equations are a fundamental building block of the mathematical description of many physical phenomena. Thus, solving problems involving complex differential equations is necessary for construction of scientific models of these phenomena, which can then be used to make useful predictions, such as weather forecasts. Aside from some simplified cases, complex differential equations cannot be solved exactly. Time integration methods are a class of numerical algorithms used to compute approximate solutions to differential equations, by stepping a given initial solution forward in time, producing a new solution at each timestep. Time integration methods are generally categorized as either explicit or implicit methods. Explicit methods are simpler, but have significant restrictions on the size of timesteps for challenging differential equations. Implicit methods relax this timestep restriction, but are much more expensive to compute. The recently introduced ``lightly-implicit'' K-methods provide a way to fuse the advantages of both implicit and explicit methods, by effectively treating a portion of the problem implicitly and the remainder explicitly. This work seeks to further two primary objectives: evaluation of these K-methods on very large problems, and development of new time integration methods. To this end, Rosenbrock--Krylov methods, the first K-methods, are applied to a large-scale problem and examined in a scalability study, and two new families of time integration methods are introduced: biorthogonal Rosenbrock--Krylov methods, and linearly implicit multistep methods. Ultimately, the goal is to develop new methods which allow for the creation of larger, more detailed, and more accurate scientific models, in order to get better and faster predictions.

Phase-field models are emerging as a promising strategy to simulate interfacial phenomena. Rather than tracking interfaces explicitly as done in sharp interface descriptions, these models use a diffuse order parameter to monitor interfaces implicitly. This implicit description, as well as solid physical and mathematical footings, allow phase-field models to overcome problems found by predecessors. Nonetheless, the method has significant drawbacks. The phase-field framework relies on the solution of high-order, nonlinear partial differential equations. Solving these equations entails a considerable computational cost, so finding efficient strategies to handle them is important. Also, standard discretization strategies can many times lead to incorrect solutions. This happens because, for numerical solutions to phase-field equations to be valid, physical conditions such as mass conservation and free energy monotonicity need to be guaranteed. In this work, we focus on the development of thermodynamically consistent algorithms for time integration of phase-field models. The first part of this thesis focuses on an energy-stable numerical strategy developed for the phase-field crystal equation. This model was put forward to model microstructure evolution. The algorithm developed conserves, guarantees energy stability and is second order accurate in time. The second part of the thesis presents two numerical schemes that generalize literature regarding energy-stable methods for conserved and non-conserved phase-field models. The time discretization strategies can conserve mass if needed, are energy-stable, and second order accurate in time. We also develop an adaptive time-stepping strategy, which can be applied to any second-order accurate scheme. This time-adaptive strategy relies on a backward approximation to give an accurate error estimator. The spatial discretization, in both parts, relies on a mixed finite element formulation and isogeometric analysis. The codes are available online and implemented in PetIGA, a high-performance isogeometric analysis framework.

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. / Doctor of Philosophy / Mathematical descriptions of physical processes are often in the form of systems of differential equations describing the time-evolution of a phenomenon. Computer simulations are realizations of these equations using well-known discretization schemes. Numerical time-stepping methods allow us to advance the state of a computer model using a sequence of time-steps. This thesis investigates current methods in time-stepping schemes, identifying a number of additional features needed to improve the speed and accuracy of simulations, and devises new methods suitable for large-scale applications where multiple processes of different physical nature drive the equations, potentially with significant differences in their time-scales. Various families of new methods are introduced with proper mathematical formulations provided for creating new ones on demand. The accuracy and stability of the methods are investigated using standard analysis techniques. These methods are then used in numerical experiments to investigate their abilities.

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. / Ph. D.

Growing interest in the simulation of first order transient systems, typical of those encountered in transient heat conduction, flow transport, and fluid dynamics, has prompted the development of a variety of time integration methods for solving these systems numerically. The primary contribution of this thesis is the design and development of a new time integration/discretization framework, under the class of single step single solve algorithms which are the most popular, for use in such first order transient systems with computationally attractive features. These include second order accuracy, unconditional stability, zero-order overshoot, and controllable numerical dissipation with a new selective control feature which overcomes the restrictions in the existing and current state-of-the-art methods. Throughout the thesis, we demonstrate the capability and advantage of the newly developed framework, termed GS4-1, in comparison to existing methods using various types of numerical examples (both linear and nonlinear). The numerical results consistently demonstrate the roles played by the new feature in improving the numerical solutions of both the primary variable and its time derivative which is important to correctly capture the dynamics of the problems, in contrast to the existing methods without such a feature. Additionally, a breakthrough contribution presented in this thesis is the development of an isochronous integration framework (iIntegrator), stemming from the novel relations between the newly developed GS4-1 framework and the existing GS4-2 framework (for second order dynamic systems). Such a development enables the use of the same computational framework to solve both first and second order dynamic systems without having to resort to the individual GS4-1 and GS4-2 frameworks; hence the practicality in the computational and implementation aspects. Finally, the application of the new GS4-1 framework to silica particle deposition, which is a practical problem of interest, is presented with the focus primarily on the physics of the problem. In this part of the thesis, a numerical model of the problem is presented and employed to investigate the effects of the flow and physicochemical parameters on the rate of deposition. The results of the parametric studies undertaken based on the employed numerical model enable some recommendations for the mitigation of the problem, and therefore serve as additional valuable contribution of the thesis.

A numerical search approach is used to design high-order diagonally implicit Runge-Kutta (DIRK) schemes equipped with embedded error estimators, some of which have identical diagonal elements (SDIRK) and explicit first stage (ESDIRK). In each of these classes, we present new A-stable schemes of orders six (the highest order of previously known A-stable DIRK-type schemes) up to order eight. For each order, we include one scheme that is only A-stable as well as one that is stiffly accurate and/or L-stable. The latter require more stages but give better results for highly stiff problems and differential-algebraic equations (DAEs). The development of the eighth-order schemes requires, in addition to imposing A-stability, finding highly accurate numerical solutions for a system of 200 equations in over 100 variables, which is accomplished via a combination of global and local optimization strategies. The accuracy, stability, and efficiency of the schemes are demonstrated on diverse problems.

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.

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. / Doctor of Philosophy / Almost all time-dependent physical phenomena can be effectively described via ordinary differential equations. This includes chemical reactions, the motion of a pendulum, the propagation of an electric signal through a circuit, and fluid dynamics. In general, it is not possible to find closed-form solutions to differential equations. Instead, time integration methods can be employed to numerically approximate the solution through an iterative procedure. Time integration methods are of great practical interest to scientific and engineering applications because computational modeling is often much cheaper and more flexible than constructing physical models for testing. Large-scale, complex systems frequently combine several coupled processes with vastly different characteristics. Consider a car where the tires spin at several hundred revolutions per minute, while the suspension has oscillatory dynamics that is orders of magnitude slower. The brake pads undergo periods of slow cooling, then sudden, rapid heating. When using a time integration scheme for such a simulation, the fastest dynamics require an expensive and small timestep that is applied globally across all aspects of the simulation. In turn, an unnecessarily large amount of work is done to resolve the slow dynamics. The goal of this dissertation is to explore new "multimethods" for solving differential equations where a single time integration method using a single, global timestep is inadequate. Multimethods combine together existing time integration schemes in a way that is better tailored to the properties of the problem while maintaining desirable accuracy and stability properties. This work seeks to overcome limitations on current multimethods, further the understanding of their stability, present new applications, and most importantly, develop methods with improved efficiency.

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. / Ph. D.

O fenômeno de impacto mecânico entre corpos sólidos está presente em diversas áreas da engenharia. Exemplos atuais deste tipo de problema podem ser encontrados no projeto de elementos de máquinas, sistemas de transporte como containers com material nuclear, tubulações em indústrias químicas, autoveículos e várias outras estruturas que devem obedecer à códigos de segurança estabelecidos por legislações governamentais. Na maioria destes casos, o conhecimento das tensões oriundas do impacto entre os corpos é fundamental para evitarem-se fa-lhas nas estruturas projetadas, predizer danos indesejáveis, diminuir coeficientes de segurança, etc. Neste contexto, é proposta neste trabalho uma metodologia de projeto contra impacto em sistemas mecânicos elásticos baseada nas equações de superfície de tensão máxima, que representam diferentes situações de impacto em uma determinada geometria. O Método dos Elementos Finitos com a integração explícita no tempo é aplicado para resolver o problema dinâmico associado ao impacto. Como exemplos de aplicações são estudados um suporte e um eixo chavetado. / Impact between solid bodies is present in many areas of engineering. Relevant examples of this sort of problem can be found in machine element design, transport systems such as containers for nuclear material, pipes in chemical plants, vehicles and many others structures that should comply with safety codes issued by govern agencies. In the majority of these cases, the knowledge of the stresses due to the impact between the bodies is fundamental to avoid failures on the designed structures, to predict undesired damages, and to decrease safety factors. Therefore, in this work a design methodology for linear mechanical systems submitted to impact is proposed. It is based on the surface of maximum stress which represents different crash situations for a given elastic model. The Finite Element Method with the explicit time integration algorithm is used to solve the associated dynamic problem. Examples are presented such as a bracket and a shaft.