Schedule for: 18w5152 - Integrating the Integrators for Nonlinear Evolution Equations: from Analysis to Numerical Methods, High-Performance-Computing and Applications

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Breakfast is served daily between 7 and 9am in the Vistas Dining Room, the top floor of the Sally Borden Building.

A brief introduction to BIRS with important logistical information, technology instruction, and opportunity for participants to ask questions.

The words "stiff", "stiffness", "stiffening", etc., arise often in applications when simulating, calibrating and controlling dynamics. But these words often have different meanings in different contexts. A subset on which we will concentrate includes: (i) Textbook-type (decaying) numerical ODE stiffness (ii) Highly oscillatory stiffness (iii) Stiffness matrix (iv) Numerical stiffening. Some of these terms are popular in scientifc computing, while others come from mechanical engineering. A potential confusion may arise in this way, and it gets serious when more than one meaning is encountered in the context of one application. Such is the case with the simulation of deformable objects in visual computing, where all of the above appear in one way or another under one roof. In this lecture I will describe the meaning of stiffness in each of these topics, how they arise, how they are related, what practical challenges they bring up, and how these challenges are handled in context. The concepts and their evolution will be demonstrated. It is about meshes -- their resolution and spectral properties -- both in time and in space.

Meanwhile, a large toolbox of numerical time integrators for nonlinear partial differential equations (PDEs) exists, based on different discretization techniques such as discretizing the variation-of-constants formula (e.g., Gautschi-type or exponential integrators) or splitting the full equation into a series of simpler subproblems (e.g., Splitting methods). In many situations these classical numerical schemes allow a precise and efficient approximation. This, however, drastically changes whenever "non-smooth" phenomena enter the scene such as for problems at low-regularity and high oscillations. Classical schemes fail to capture the oscillatory parts within the solution which leads to severe instabilities and loss of convergence. In this talk I present a new class of Fourier-integrators for dispersive equations at low-regularity and high oscillations. The key idea in the construction of the new schemes is to tackle and deeply embed the underlying structure of resonances in the numerical discretizations. These are the cornerstones of theoretical analysis of the long time behaviour of differential equations and their numerical discretizations (cf. Modulated Fourier Expansion; Hairer, Lubich & Wanner) and will offer the new schemes strong geometric structure at low regularity.

State-of-the-art research within the field of Visual Computing enables the acquisition, creation, processing, and manipulation of visual content like images and 3D geometry. The investigated computational methods and algorithms allow and facilitate applications beyond classical computer graphics and animation, such as computer aided product design and fabrication as well as computational architecture, the creation of synthetic training data for machine learning, special effects in movies and interactive applications hardly distinguishable from reality, tomography and medical imaging, and computational photography. This often requires the appropriate handling of the underlying physics focusing mostly on globally accurate simulations by providing numerical tools that intrinsically respect key defining properties of the physical systems. At the same time, large complexities, the accurate coupling and interaction of different types of physical systems as well as the interactions within these systems have to be addressed properly. Moreover, several applications require interactivity defining hard constraints with respect to acceptable computation times. This talk provides an overview of different aspects of Visual Computing, and presents a selection of the speaker's work within this field as concrete examples in order to illustrate potential overlap of this application-oriented discipline with the numerical mathematics community aiming for the stimulation of discussions about mutual research interests and potential collaborative future work.

The efficient numerical simulation of stiff multiphysics systems remains a core challenge in scientific computing. Examples are fluid structure interaction, earth system models or turbulent flames. We consider problems with the following characteristics: They are large scale, all components are stiff, possibly on different time scales and there are codes for the subproblems available. Thus, we want a partitioned numerical method, meaning that reuse of the existing codes is possible. Thereby, we assume that while we have access to the source codes, we want to edit that code as little as possible. In particular we assume that we can repeat a time step. We are then looking for numerical methods that are implicit and at least order two, time adaptive, allow the subsolvers to run in parallel and allow for different time steps in the different models. We are not aware of a method that fulfills all of these properties and suggest two methods of our own for the case of two systems being coupled. The core idea is the following: We have a time integration method of at least order two for each subproblem and assume that we can restart these with new initial data and that during time integration, information for the other solver at all times can be provided using interpolation. This continuous representation of the numerical solution is updated after each local time step. Then the solvers run in parallel over a macro time window and are free to choose their own timesteps in an adaptive way without outside interference. At the end of the macrostep, it is checked if the coupled system is fulfilled up to a tolerance, if not, the time window is repeated. Crucial questions are order of the time integration method and convergence of the time window iteration, also called waveform relaxation. This is shown numerically for representative test cases. For the specific case of two linear heat equations with different material properties coupled across an interface, we suggest to do the waveform relaxation in the form of a Neumann-Neumann coupling, known from domain decomposition. There, the choice of the relaxation parameter is crucial and previous analysis by Gander and Kwok for the semidiscrete case does not apply. We thus perform a fully discrete analysis for the case of fixed but different time steps for the subproblems. Numerical results show that this can be used for the time adaptive case as well.

Lunch is served daily between 11:30am and 1:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Meet in the Corbett Hall Lounge for a guided tour of The Banff Centre campus.

Meet in foyer of TCPL to participate in the BIRS group photo. The photograph will be taken outdoors, so dress appropriately for the weather. Please don't be late, or you might not be in the official group photo!

In this talk we discuss a Bayesian approach for inverse problems involving elliptic differential equations with multiple scales. Computing repeated forward problems in a multiscale context is computationnally too expensive and we propose a new strategy based on the use of "effective" forward models originating from homogenization theory. Convergence of the true posterior distribution for the parameters of interest towards the homogenized posterior is established via G-convergence for the Hellinger metric. A computational approach based on numerical homogenization and reduced basis methods is proposed for an efficient evaluation of the forward model in a Markov Chain Monte Carlo procedure. We also discuss a methodology to account for the modeling error introduced by the effective forward model and the combination of the Bayesian multiscale method with a probabilisitic approach to quantify the uncertainty in building the effective forward model for a multiscale elastic problem in random media. References: A. Abdulle, A. Di Blasio, Numerical homogenization and model order reduction for multiscale inverse problems, to appear in SIAM MMS. A. Abdulle, A. Di Blasio, A Bayesian numerical homogenization method for elliptic multiscale inverse problems, Preprint submitted for publication.

Efficient software packages are often the main vehicle for inserting numerical methods developed in the applied mathematics community into complex scientific simulation codes. However, strong assumptions on use contexts or computing system architectures in the user interfaces can make the difference between a package providing useful methods or providing reasons to reject input from the mathematics computing. Flexible packages with effective user interfaces greatly ease the transition of new methods into scientific software. SUNDIALS is a suite of robust and scalable solvers for systems of ordinary differential equations, differential-algebraic equations, and nonlinear equations designed for use on computing systems ranging from desktop machines to super computers. The suite consists of six packages: CVODE(S), ARKode, IDA(S), and KINSOL, each built on common vector and solver application programming interfaces (API) allowing for application-specific and user-defined linear solvers, nonlinear solvers, data structures, encapsulated parallelism, and algorithmic flexibility. In this presentation we will overview the design principles adopted by the SUNDIALS development team and discuss how they are manifested in package flexibility and user interfaces. In addition, we will overview the current interfaces in SUNDIALS with examples demonstrating how the interfaces work in applications and the benefits to both package developers and users. This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. LLNL-ABS-760560.

In this lecture, we review the existing and emerging time integration practices used in the operational NWP models. We will emphasize the reasons why such numerical strategies were adopted and why many others have been disregarded.

In this session graduate students and postdocs participating in the workshop briefly outline their research.

A buffet dinner is served daily between 5:30pm and 7:30pm in the Vistas Dining Room, the top floor of the Sally Borden Building.

Computers will soon reach exascale performance, i.e. the capability to perform more than 10^18 FLOPS. This will e.g. enable flow simulations where many suspended particles can be resolved as geometric objects and the when time evolution is computed coupling the hydrodynamics of the fluid phase with rigid body dynamics. Or, on a larger physical scale, the whole volume of the Earth can be discretized with a global mesh of 1km resolution to simulate Earth mantle dynamics. However, such extreme scale simulations will only be possible with parallel algorithms that can exploit millions of processor cores. In this talk we will report early experience with such massively parallel computations for multiphysics problems.

I will discuss two topics: characteristic-based flux partitioning and a posteriori error estimation. In the first part I introduce a characteristic-based flux partitioning for the semi-implicit time integration of atmospheric flows and discuss it in the context of compressible Euler equations. Here the acoustic time-scale is significantly faster than the advective scale, yet it is typically not relevant to atmospheric and weather phenomena. The acoustic and advective components of the hyperbolic flux are separated in the characteristic space. Results with high-order, conservative additive Runge-Kutta methods are briefly discussed. In the second part I will briefly discuss some new time-stepping strategies with built-in global error estimators. These methods can be cast as general linear schemes that provide pointwise a posteriori errors. I will show some preliminary results on ODE and PDE problems.

We consider additively partitionable time dependent problems of the form $y'(t) = f_1(y) + f_2(y) + ...$, where there is considerable difference in the underlying time scales of each operator. Many important application problems are of this form, such as reaction-transport models for combustion and stellar dynamics, climate simulation, etc. To minimize computational work, it is desirable to time step the overall problem in a “multirate” manner using different time step sizes for each operator while still preserving overall accuracy and stability. The order of accuracy, stability, and overall computational efficiency of a multirate method depend intimately on how the coupling between partitions is constructed. The coupling structure is affected by the ratio of fast to slow step sizes and has a high number of degrees of freedom when the fast-to-slow ratio is large. In this talk, we consider some of the resulting difficulties in developing stable implicit multirate methods and discuss some initial progress in constructing effective schemes.

Although originally constructed as a library for additive Runge-Kutta methods, we have recently overhauled the ARKode library to serve as infrastructure for general, adaptive, one-step time integration methods. Sharing a common set of vector, matrix, linear solver and nonlinear solver objects, this new infrastructure enables the rapid development of novel time integration algorithms for the parallel solution of large-scale ODE systems. In this talk, I will discuss this updated infrastructure, as well as our new results in constructing solvers for methods of both GARK (ImEx or fully explicit) and MIS (multirate with explicit slow component) type.

The idea of preconditioning iterative methods for the solution of linear systems goes back to Jacobi (1845), who used rotations to obtain a system with more diagonal dominance, before he applied what is now called Jacobi's method. The preconditioning of linear systems for their solution by Krylov methods has become a major field of research over the past decades, and there are two main approaches for constructing preconditioners: either one has very good intuition and can propose directly a preconditioner which leads to a favorable spectrum of the preconditioned system, or one uses the splitting matrix of an effective stationary iterative method like multigrid or domain decomposition as the preconditioner. Much less is known about the preconditioning of non-linear systems of equations. The standard iterative solver in that case is Newton's method (1671) or a variant thereof, but what would it mean to precondition the non-linear problem ? An important contribution in this field is ASPIN (Additive Schwarz Preconditioned Inexact Newton) by Cai and Keyes (2002), where the authors use their intuition about domain decomposition methods to propose a transformation of the non-linear equations before solving them by an inexact Newton method. Using the relation between stationary iterative methods and preconditioning for linear systems, we show in this presentation how one can systematically obtain a non-linear preconditioner from classical fixed point iterations, and present as an example a new two level non-linear preconditioner called RASPEN (Restricted Additive Schwarz Preconditioned Exact Newton) with substantially improved convergence properties compared to ASPIN.

Differential equations arising in many practical applications are characterized by multiple time scales. Multirate time integration seeks to solve them efficiently by discretizing each scale with a different, appropriate time step, while ensuring the overall accuracy and stability of the numerical solution. While the multirate idea is elegant and has been around for decades, multirate methods are not yet widely used in applications. This is due, in part, to the difficulties raised by the construction of high-order multirate schemes. We discuss the design of practical high-order multirate methods using the theoretical framework of generalized additive Runge--Kutta methods MR-GARK, which provides the generic order conditions and the linear and nonlinear stability analyses. In a seminal paper Knoth and Wolke (1998) proposed a hybrid solution approach: discretize the slow component with an explicit Runge-Kutta method, and advance the fast component via a modified fast differential equation. The idea led to the development of multirate infinitesimal step (MIS) methods in (Wensch et al. 2009). Guenther and Sandu (2016) explained MIS schemes as a particular case of MR-GARK methods. The hybrid approach offers extreme flexibility in the choice of the numerical solution process for the fast component. This work discusses new families of multirate infinitesimal GARK schemes (MRI-GARK) that extends the MIS approach in multiple ways. Order conditions theory and stability analyses are developed, and practical explicit and implicit methods of up to order four are constructed. Numerical results confirm the theoretical findings. We expect the new MRI-GARK family to be most useful for systems of equations with widely disparate time scales, where the fast process is dispersive, and where the coupling between the fast and slow dynamics is relatively weak.

Krylov subspace spectral (KSS) methods are high-order accurate, explicit time-stepping methods with stability characteristic of implicit methods. This "best-of-both-worlds" compromise is achieved by computing each Fourier coefficient of the solution using an individualized approximation, based on techniques from "matrices, moments and quadrature" due to Golub and Meurant for computing bilinear forms involving matrix functions. The result is superior scalability to that of other time-stepping approaches, which motivates continued development of KSS methods for high-resolution simulation. Through combination with EPI methods due to Tokman, et al., KSS methods have been shown to be applicable to nonlinear PDEs as well. This talk will present an overview of their derivation and essential properties, including new theoretical results, and also highlight ongoing projects aimed at enhancing their performance and applicability.

In this talk I will describe and demonstrate two software packages I’ve developed as tools for research on time integration methods: Nodepy and RK-Opt. RK-Opt is a MATLAB package for optimizing properties of Runge-Kutta methods. It uses nonlinear optimization to find methods with optimized accuracy and stability properties with a given order, number of stages, and structure (explicit, implicit, diagonally implicit, low storage, etc.). It can be used to design specialized integrators tailored specifically to a class of PDE semi-discretizations. Nodepy (Numerical ODEs in Python) is a package for studying the properties of time integrators, including Runge-Kutta, linear multistep, extrapolation, and deferred correction methods. Given the coefficients of a method, it can provide a wealth of information on its accuracy and stability properties. It also includes basic tools for investigating a method’s efficiency in practical terms.

I will give a brief overview of iterative methods based on Spectral Deferred Corrections and highlight how the versatility of the SDC approach allows for the straightforward construction of specialized integration schemes for complex physical problems with high-order of accuracy. I will then explain how SDC methods have been incorporated into the parallel-in-time PFASST algorithm and highlight some recent successes.

The current trend in high-performance computing of a stagnating or even decreasing processor speed poses new challenges to solve PDEs within a particular time frame. Here, disruptive mathematical reformulations which exploit additional degrees of parallelism also in the time dimension gained increasing interest over the last two decades. One of such reformulations is a rational approximation of exponential integrators (REXI) which replaces a purely CFL-limited and therefore sequential time integration for linear oscillatory or diffusive systems by a sum of solutions of decoupled systems. Each of the terms in the sum can then be solved independently, hence massively parallel for arbitrarily long time step sizes of linear operators. We will present studies conducted with the linear and non-linear shallow-water equations on the rotating sphere. These equations represent a simplified model of the equations used for weather and climate simulations, but with a focus on horizontal aspects. Various benchmarks will be discussed, including timestepsize-to-error and wallclocktime-to-error, revealing sweet spots of exponential integrators. The results motivate further explorations of REXI for operational weather/climate systems.

In this work we present recent developments of multiderivative time integrators for time dependent ordinary and partial differential equations. Such methods are an umbrella class of solvers that contain all Taylor (Lax-Wendroff) as well as Runge-Kutta (multistage) solvers. These methods increase the orders of accuracy of the base solvers via two mechanisms: 1) adding stage values, and 2) incorporating higher derivatives of the unknown. We include results that leverage differential transforms to define higher derivatives of the unknown variables, embedding these solvers within the spectral deferred correction framework, as well as definitions for strong stability preserving versions of these solvers.

Adaptive mesh refinement (AMR) has been widely used for solving a variety of PDEs in which solution features of interest are spatially localized. In the past three decades, much attention has been devoted to spatial discretization in the AMR setting, while relatively little attention has been paid to time stepping on a multi-resolution hierarchy of grids. Adaptive time stepping, aimed at preserving a constant CFL number across the mesh hierarchy, was built into the original Berger-Oliger and Berger-Collela approaches to adaptive mesh refinement, but in its original form, this adaptive time stepping was only well suited to single-step, single-stage time stepping schemes. While more sophisticated schemes can be easily implemented if the same size time step is used on all grids (global time stepping), little attention has been given to combining adaptive time stepping with multi-step or multi-stage schemes. In this talk, we present our current efforts to combine the multi-stage, explicit Runge-Kutta-Chebyshev (RKC) integration scheme with adaptive time stepping. Using the RKC scheme, we can solve parabolic equations using explicit time stepping in a manner that is competitive, if not superior to alternative implicit approaches. The multi-stage nature of the scheme requires special algorithmic structures be built into the core AMR time stepping strategy. In particular, for finite volume schemes, extra effort may be required to ensure that the resulting scheme is conservative. We have developed a algorithmic approach to implementing the RKC scheme on an one-dimensional adaptively refined mesh using adaptive time stepping, and are currently working on a two dimensional approach in our block-structured adaptive code ForestClaw. We will present results from typical reaction diffusion problems.

Heterogeneous multiscale algorithms are methods designed to cope with the challenge of solving numerically problems with different temporal or spatial scales. The talk considers the case of delay ordinary differential equations subject to periodic forcing in cases where the solution has to be computed over long time intervals. We shall show two alternative techniques; one of them may achieve arbitrarily high orders of convergence with a computational complexity that is independent of the forcing frequency.

In the numerical solution of the algebraic Riccati equation $A^* X + XA − XBB^∗ X + C^∗ C = 0$, where $A$ is large, sparse and stable, and $B$, $C$ have low rank, projection methods have recently emerged as a possible alternative to the more established Newton-Kleinman iteration. A robust implementation of these methods opens to new questions on the use of dissipativity properties of the given matrix $A$. In this talk we briefly discuss the algorithmic aspects of projection methods, together with some new hypotheses that ensure their well posedness. If time allows, considerations on the differential Riccati equation will be included.

We study the unconditional stability of implicit–explicit BDF methods applied to Dahlquist’s first test equation. This can be viewed as an attempt to extend the well–known $A(\theta)$-stability concept to implicit–explicit multistep methods.

We present the construction and analysis of uniformly accurate nested Picard iterative integrators (NPI) for the Dirac equation in the nonrelativistic limit involving a dimensionless parameter inversely proportional to the speed of light. To overcome the difficulty induced by the rapid temporal oscillation, we present the construction of several NPI methods which are uniformly first-, second- and third-order convergent in time. The NPI method can be extended to arbitrary higher order in time with optimal and uniform accuracy. The implementation of the second order NPI method will be demonstrated and analyzed.

Many problems encountered in plasma physics require a kinetic description. The associated partial differential equations are posed in an up to six-dimensional phase space. A direct discretization of this phase space, often called the Eulerian approach, has many advantages but is extremely expensive from a computational point of view. In this talk we propose a dynamical low-rank approximation to the Vlasov equation. This approximation is derived by constraining the dynamics to a manifold of low-rank functions via a tangent space projection. Then the projection is split into the sub-projections from which it is built. This reduces a time step for the six- (or four-) dimensional Vlasov--Poisson equation to solving two systems of three- (or two-) dimensional advection equations. This projector-splitting approach also enables us to dynamically adjust the rank during the simulation.

In this talk, we will describe some fractional step (time-spitting) strategies for the Navier-Stokes, Stokes and stationary Stokes coupled to transport-diffusion-reaction. The main objective is that the resulting PDEs can be solved numerically by fast solvers (particle methods, fast interpolation, FFT). It turns out that the flow computation itself doesn't require sophisticated numerical analysis since the data are structured, but the splitting of transport-reaction-diffusion bumps into some limitations due to the time multi-scale feature and stiffness. One of the most interesting perspectives is the use of exponential integrators in order to handle different operators at the same time and blend it together with particle-strength-exchange schemes.

In this talk, we consider the numerical solution of the highly-oscillatory Vlasov equations. Designed in the spirit of recent uniformly accurate methods, the scheme remains insensitive to the stiffness of the problem in terms of both accuracy and computational cost. The method is based on a careful ad-hoc reformulation of the equations. Some numerical results will be given to illustrate the behavior of the method.

In this talk, we will present some new results on the backward error analysis for the approximation of the action of the matrix exponential $\exp(A)v$. Al-Mohy and Higham (Computing the action of the matrix exponential with an application to exponential integrators, SIAM J. Sci. Comput. 33:488-511, 2011) endowed the scaled Taylor method with a backward error analysis approach based on the norm of the matrix $A$ and of its first few powers. Such an approach allows to choose the scaling parameter $s$ and the polynomial degree of approximation $m$ granting the desired accuracy. C., Kandolf, Ostermann, and Rainer (The Leja method revisited: backward error analysis for the matrix exponential, SIAM J. Sci. Comput. 38:A1639-A1661, 2016) extended it to Leja points interpolation. Moreover, they considered a new contour integral expansion of the backward error onto an ellipse $\Gamma=\partial K$, such that $$ \Lambda_{\varepsilon}(A)\subseteq W(A)+\Delta_\varepsilon\subseteq K, $$ where $\Lambda_\varepsilon(A)$ is the $\varepsilon$-pseudo-spectrum of $A$, $W(A)$ is its field of values and $\Delta_\varepsilon$ is the disk of radius $\varepsilon$, with a \emph{fixed} value of $\varepsilon$. Such an approach is particularly suited for matrices with skinny pseudo-spectra, such as diffusion, transport, and Schrödinger matrices, since it better adapts to the matrix and mitigates phenomena like under-scaling or over-scaling, leading to an improved accuracy and a lower computational cost. Therefore, we decided to enhance this technique by optimizing the choice of $\varepsilon$ (which now depends on the matrix $A$) and by using the so-called Leja-Hermite points introduced by C., Kandolf and Zivcovich (Backward error analysis of polynomial approximations for computing the action of the matrix exponential, BIT Numer. Math., 2018). This is a joint work with Franco Zivcovich.

Multirate time integration methods allow for adaptation of local time step size, are thus particularly efficient for multi-scale systems that have widely different time scale characteristics. However, classic multi-rate schemes designed for ODEs are typically difficult to be extended for PDEs because of additional stability requirements. For example, hyperbolic conservation laws require the numerical schemes to be positivity-preserving and total variation diminishing. While multi-rate partitioned Runge-Kutta methods were proposed to use different time steps on different grid cells to overcome the global CFL constraint, buffer regions have to be introduced in order to accommodate the transition between subdomains. Implementation of these methods are tricky and highly application-dependent, and especially challenging for general-purpose libraries. In this talk, we will present our recent development of multi-rate time integrators in PETSc, which is one of the most popular software library for the scalable solution of scientific applications modeled by PDEs. We will introduce the new data structure support for componentwise partitioned systems that arise naturally from mesh refinement, and show how the serial and parallel efficiency is achieved. To illustrate the capability of the PETSc multi-rate integrators for hyperbolic conservation laws, we will focus on the advection equation solved using conservative finite volume methods with consistent slope limiters.

A proof of convergence is given for semi- and full discretizations of mean curvature flow of closed two-dimensional surfaces. The numerical method proposed and studied here combines evolving finite elements, whose nodes determine the discrete surface like in Dziuk's method, and linearly implicit backward difference formulae for time integration. The proposed method differs from Dziuk's approach in that it discretizes Huisken's evolution equations for the normal vector and mean curvature and uses these evolving geometric quantities in the velocity law projected to the finite element space. This numerical method admits a convergence analysis, which combines stability estimates and consistency estimates to yield optimal-order $H^1$-norm error bounds for the computed surface position, velocity, normal vector and mean curvature. The stability analysis is based on the matrix--vector formulation of the finite element method and does not use geometric arguments. The geometry enters only into the consistency estimates. Numerical experiments illustrate and complement the theoretical results.

5-day workshop participants are welcome to use BIRS facilities (BIRS Coffee Lounge, TCPL and Reading Room) until 3 pm on Friday, although participants are still required to checkout of the guest rooms by 12 noon.