Schedule for: 16w5030 - Fifth Parallel-in-time Integration Workshop

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.

Domain decomposition methods like the classical Schwarz method and the Dirichlet-Neumann and Neumann-Neumann methods have historically been developed for steady partial differential equations. All these methods have however also a natural waveform relaxation extension to time dependent partial differential equations. These waveform relaxation variants are still based on a spatial decomposition of the physical domain into subdomains, but then time dependent problems are solved in the subdomains, and information is exchanged on the space-time interfaces between subdomains in an iteration that converges in space-time to the underlying solution of the evolution problem. While domain decomposition methods for steady problems converge linearly, the waveform relaxation variants exhibit superlinear convergence for parabolic problems, and often become direct solvers for hyperbolic problems. After a simple introduction to the main classes of domain decompositions for steady problems, I will give in my talk an overview over the development of their waveform relaxation variants as it happened over the past 20 years.

While it has been established that Parareal has stability problems for hyperbolic and advection-dominated problems, details of how it propagates waves are less well understood. The talk will show how, starting from the interpretation of Parareal as a preconditioned fixed point iteration, one can derive a stability function for linear problems. After ``normalising'' this function to a unit time interval it is then possible to derive and analyse a discrete dispersion relation for Parareal. This allows to estimate the impact of e.g. the choice of propagators, size of fine and coarse time step, number of time slices etc. on Parareal's wave propagation characteristics. In particular, I will discuss the effects of phase and amplitude errors in the coarse method on convergence. The formulation also allows a worst case estimate for convergence through the maximum singular value of the error propagation matrix. This allows to link the number of iterations to the number of processors in the speedup model and to make better predictions about weak scaling of Parareal.

For time-dependent PDEs, parallel-in-time integration using the "parallel full approximation scheme in space and time" (PFASST) is a promising way to accelerate existing space-parallel approaches beyond their scaling limits. While many use cases and benchmarks exist, a solid and reliable mathematical foundation is still missing. In this talk, we formulate PFASST as a specialized FAS multigrid method. We use spectral deferred corrections for the definition of block smoothers and define the appropriate coarse grid correction to establish a tight link between PFASST and standard multigrid methods, providing an easy access to the mathematical analysis and algorithmic optimization. Using local Fourier analysis, we describe first steps towards a semi-algebraic convergence analysis for the linear case and show some results for diffusive and advective prototype problems.

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. Please don't be late, or you will not be in the official group photo! The photograph will be taken outdoors so a jacket might be required.

The need for parallel-in-time is being driven by changes in computer architectures, where future speedups will be available through greater concurrency, but not faster clock speeds, which are stagnant. This leads to a bottleneck for sequential time marching schemes, because they lack parallelism in the time dimension. In this talk, we examine an optimal-scaling parallel time integration method, multigrid reduction in time (MGRIT). MGRIT applies multigrid reduction to the time dimension by solving the (non)linear systems that arise when solving for multiple time steps simultaneously. The result is a versatile approach that is nonintrusive and wraps existing time evolution codes. In this talk, we present the MGRIT framework and then discuss some recent theoretical convergence results. Some typical simplifying assumptions are made, such as a two-grid method, uniform time-line and constant coefficient PDE. The convergence estimates are then explored in a variety of settings, including common parabolic and hyperbolic spatial discretizations coupled with implicit Runge-Kutta time-stepping schemes. Overall, the convergence estimates are sharp when compared to numerical experiments and are good predictors of multilevel results. LLNL-ABS-695327

RIDC software workshop

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.

The "parallel full approximation scheme in space and time" (PFASST) is an iterative, multilevel strategy for the temporal parallelization of ODEs and discretized PDEs. In numerous studies, this approach has been successfully coupled to space-parallel solvers which use finite differences, spectral methods or even particles for discretization in space. In this talk, we will report on our experience using PFASST in time together with finite elements in space. In particular, we discuss modifications necessary to treat the mass matrix appropriately on all levels of the space hierarchy and describe a procedure to switch between these levels. For our experiments, the base implementation is the PFASST++ software, a C++ implementation of PFASST and its building blocks MLSDC and SDC. In the context of the finite element discretizations we make use of the "Distributed and Unified Numerics Environment" (DUNE), which is a modular framework for solving PDEs with grid-based methods. Using this coupling of PFASST++ and DUNE, we will show first results for a nonlinear two-component Gray-Scott model. This work is conducted within the project "ParaPhase", where the final goal is the development of a highly scalable space-time parallel adaptive algorithm for the simulation of phase-field models.

The stagnated increase in CPU frequency and the resulting trend in HPC towards massively parallel systems poses new challenges for HPC applications. The ongoing trend towards massive parallelism (accelerator cards and many-core systems) requires redesigning existing algorithms to cope with these architectures. In this presentation we will focus on oscillatory problems. Solving them with a regular time stepping method leads to a timestep-by-step sequentialization in the time dimension. In combination with the CFL limitation, an increasing resolutions also results in an increase in the number of these sequential time steps. For problems which are already limited in their scalability this directly leads to an increase in wall clock time which can make the results less valuable or even useless. Based on underlying properties of oscillatory problems we make use of a recently developed method of a "rational approximation of exponential integrator" (REXI). This method uses features of linear oscillatory problems which allows a massively parallel formulation. This yields several beneficial advantages: e.g. an increase in resolution does not lead to an increase in the number of time steps as it is the case with conventional time stepping methods. In contrast, it leads to an increase in the degree of parallelization. Our results show speedups of over 2 orders of magnitude compared to a standard time stepping method and a scalability model indicates robust scalability beyond 100k cores in case of large time step sizes. Finally, an extension to a semi-Lagrangian time stepping scheme to cope with non-linearities will be discussed.

Unsteady turbulent flow simulations using the Navier Stokes equations require larger and larger problem sizes. On an other side, new supercomputer architectures will be available in the next decade, with computational power based on a larger number of cores rather than significantly increased CPU frequency. Hence most of the current generation CFD software will face critical efficiency issues if bounded to massive spatial parallelization and we consider time parallelization as an attractive alternative to enhance efficiency on multi- cores architectures. Several algorithms developed in the last decades (Parareal, PFASST) may be straightforwardly applied to the Navier-Stokes equations, but the Parareal algorithm remains one of the simplest solutions in the case of explicit time stepping, compressible flow Based on an optimized implementation of Parareal,3 we modelize the speed-up obtained when combining both space and time parallelizations. This modelization takes into account the speedup of an actual structured, massively parallel CFD solver and the cost of time communications, both measured on two different supercomputers. Some preliminary requirements for a worthy time-parallel integration will be then derived, in terms of both Parareal iteration count and size of the time subdomain window. We then study within this framework, possible enhancements of the well-known convergence difficulties for Parareal encountered for advection dominated problems. The proposed approach is based on the representation of Parareal as an algebraic system of nonlinear equations solved by a preconditioned Newton’s method. The new formulation targets the reduction of the degree of non-normality of its Jacobian by slightly modifying the Parareal iteration. Performance on examples related to canonical linear problems, like the Dahlquist and the one- dimensional advection equation, is analysed. To conclude we comment on the extension of this method to nonlinear problems.

Participants are invited to pose open questions and/or topics for discussion to the rest of the group. Those wishing to present questions should prepare quick presentations (2-3 slides) for context. Participants are encouraged to discuss these questions openly throughout the remainder of the workshop.

In recent years we have developed reliable order $p$ methods for the approximate solution of general systems of IVPs and DDEs. We have used these methods to implement effective techniques for parameter estimation and sensitivity analysis for IVPs and DDEs. The two techniques we have investigated are the "variational approach" and the "adjoint approach". We will discuss the cost and reliability of each approach and identify several factors that contribute to the performance of each. In particular we will discuss why, for each technique, the underlying numerical IVP or DDE solver must exploit inherent parallelism if the number of state variables, the number of parameters, or the number of data points is large.

In recent years we have developed reliable order $p$ methods for the approximate solution of general systems of IVPs and DDEs. We have used these methods to implement effective techniques for parameter estimation and sensitivity analysis for IVPs and DDEs. The two techniques we have investigated are the "variational approach" and the "adjoint approach". We will discuss the cost and reliability of each approach and identify several factors that contribute to the performance of each. In particular we will discuss why, for each technique, the underlying numerical IVP or DDE solver must exploit inherent parallelism if the number of state variables, the number of parameters, or the number of data points is large.

We discuss an application of PFASST to simulate climate on a sphere. For the shallow water equations, the High-Order Methods Modeling Environment (HOMME) with spectral elements is used. Time parallelism on the other hand is realized by means of LIBPFASST. The focus is on implementation as well as efficiency.

Solving industrial sized problems in the field of computational fluid dynamics is a big challenge. Especially in flows with high Reynolds numbers (flows with a high ratio between convection and diffusion) small turbulent structures have to be resolved with a very fine grid in space and time. These fine grids lead, even with spatial parallelization, to long computational times, since also a reasonable amount of physical time has to be simulated. This problem can be partially circumvented by modelling the turbulent structures as it is done with statistical turbulence models. These modelling techniques allow coarser grids and shorter simulation runtimes with the drawback of a less accurate solution. Nevertheless, these simulations can still have a runtime in the range of months. For these problems a parallelization in time in addition to the parallelization in space is very appealing. The parallelization in time can be used to reduce the computational load of the typically longer physical time which is to be simulated. The easiest parallel in time method, the Parareal method, would be a good starting point for the runtime reduction. Unfortunately, it was already shown by multiple publications, that the Parareal method in its original formulation is not suitable for high Reynolds flows [e.g. Kreienbuehl et al., 2015]. The Parareal method has stability problems which occur with high Reynolds number flows. These stability problems can be investigated by applying the Parareal method to the viscous Burgers equation, which is closely related to the Navier-Stokes equations. Therefore, a parameter study was done by varying the Reynolds number of the simulated flow. In addition, it was investigated whether a scale difference of the structures in the flow have an impact on the convergence. These structures which model the turbulent scales are induced according to [Kooij et al., 2015]. Furthermore, the effect of applying different time stepping schemes was studied. The used schemes were of explicit Runge-Kutta and IMEX Runge-Kutta nature. Moreover, it will be presented whether it is possible to reduce the stability problems of applying the Parareal method to high Reynolds number flows with a semi-Lagragian formulation of the examined equations. The work of Andreas Schmitt is supported by the 'Excellence Initiative' of the German Federal and State Governments and the Graduate School of Computational Engineering at Technische Universit\"at Darmstadt.

In this paper we present an adjoint solver for the multigrid in time software library XBraid. XBraid provides a non-intrusive approach for simulating unsteady dynamics on multiple processors while parallelizing not only in space but also in the time domain. It applies an iterative multigrid reduction in time algorithm to existing spatially parallel classical time propagators and computes the unsteady solution parallel in time. However, in many engineering applications not only the primal unsteady flow computation is of interest but also the ability to compute sensitivities that determine the influence of design changes to some output quantity. In recent years, adjoint solvers have widely been developed which propagate sensitivity information backwards through the time domain. We develop an adjoint solver for XBraid that enhances the primal iterations by an iteration for computing adjoint sensitivities. In each iteration, the adjoint code runs backwards through the primal XBraid actions and computes the consistent discrete adjoint sensitivities parallel in time. It is highly non-intrusive as existing adjoint time propagators can easily be integrated through the adjoint interface. We validate the adjoint code by applying it to an unsteady partial differential equation that mimics the behavior of separated flows past bluff bodies. In our 1D model, the near wake is mimicked by a nonlinear ODE, namely the Lorenz attractor which exhibits self-excited oscillations. The far wake is modeled by an advection - diffusion equation whose upstream boundary condition is determined by the ODE mimicking the near wake. We demonstrate the integration of a serial time stepping algorithm, that solves the PDE forward in time, into the parallel-in-time XBraid framework as well as the development of the corresponding adjoint interface. The resulting sensitivities are in good agreement with those computed from finite differences. Nevertheless, there is still great potential for optimizing the performance of the adjoint code using advanced algorithmic differentiation techniques such as reverse accumulation and checkpointing. Due to the iterative nature of the primal and the adjoint flow computation, the method is very well suited for simultaneous optimization algorithms like the One-shot method which solve the optimization problem in the full space. They have proven to be very efficient for optimization with steady-state PDE constraints while its application to unsteady PDE is still under development. The non-intrusive adjoint XBraid solver is therefore highly desirable and will extend its application range from pure simulation to optimization.

When solving optimal control problems over a long time horizon, one can introduce additional parallelism in time by subdividing the time horizon into smaller, non-overlapping time intervals and by solving these subproblems in parallel. If the intermediate state and adjoint between time intervals are known exactly, this procedure yields the exact solution. Thus, the problem reduces to solving a nonlinear system in these intermediate states, which are related via certain propagation operators. In this talk, we present a parareal approach for solving this nonlinear system: here, the global problem is approximated by a simpler one using coarse propagators, while the fine propagation is performed in parallel over different time intervals. One then iterates until the intermediate states are consistent across time intervals. Unlike parareal for initial value problems, the coarse problem still contains a forward-backward coupling, but it is much cheaper to solve than the global fine problem. We analyze the convergence of the new method for a model linear problem and illustrate its behaviour numerically for nonlinear problems in which the control enters as an additive source term.

Recently several waveform relaxation algorithms have been proposed, where the standard Dirichlet or Robin coupling conditions are replaced by transmission conditions that are implemented in stages, for example the Dirichlet--Neumann coupling conditions. These new transmission conditions promise better rates of convergence at the expense, as it turns out, of parallel efficiency. This talk focuses on a practical implementation of these new waveform relaxation algorithms, as well as an analysis of the communication cost and efficiency overhead of these algorithms.

We consider parabolic partial differential equations defined on multiple domains. These domains are coupled at the boundary by Robin boundary conditions where the region of overlap is time-dependent. The rate of change of the geometry is much faster than the time scale of heat conduction. We apply the spectral deferred correction approach as well as a splitting in slow and fast components to this type of problem. We use the PFASST concept to obtain a parallel implementation of these concepts. One basic ingredient of PFASST are the underlying spectral deferred correction methods. Spectral deferred correction (SDC) methods start from a provisional solution. Using a simpler time integrator, this provisional solution will be iteratively improved. There are several adaptions of SDC for multi-rate problems. MISDC methods treat every scale independently in the sweeper and allows to construct high-order multi-rate methods. Typically, slow processes are treated explicitly and fast processes implicitly. We adapt the idea of the MISDC methods for the coupled heat equation and treat the diffusion part implicitly, but the fast sources in an explicit manner to avoid implicit solves for the geometry. The method will be analyzed with respect to order and stability. Finally, we present numerical results.

In this talk we discuss stable space-time discretization schemes, which allow the use of standard parallel finite element libraries to solve the arising space-time linear systems efficiently. Moreover we will have a look at space-time multigrid methods which also allow parallelization with respect to space and time.

In this study, we introduce the reduced basis methods (RBMs) to improve a convergence rate of the Parareal method for hyperbolic PDEs. We extract a small subspace consisting of main modes that compose the accurate solution from the data calculated by the fine solver during iterations. Once we got a set of reduced basis, the computational cost of time marching becomes low because of the small subspace, and therefore we can use a fine time step width. Then we can perform a coarse solver with the time step width same as a fine solver. Thus, it is expected that the convergence can be improved and the computation cost can be reduced. However, if the subspace cannot be reduced to small one, the RBMs maybe does not work well, e.g., for the complex phenomena such as the turbulence flow. We need to know what case the RBMs work well for, but currently it is not clear. Therefore, we investigate the convergence characteristics of the Parareal method for hyperbolic PDEs which have different complexity with the RBMs. In this presentation, we discuss about the results of convergence of the Parareal method with the RBMs for the linear advection equation, viscous Burgers' equation and Navier-Stokes equations.

Since clock speeds are no longer increasing, time integration is becoming a sequential bottleneck. The multigrid reduction in time (MGRIT) algorithm is an approach for exploiting parallelism in the time dimension that is designed to build on existing codes and time integration techniques. The XBraid library is an open source implementation of MGRIT. One important technique used by current simulation codes is adaptivity in both space and time. In this talk, we discuss approaches taken in XBraid to support adaptivity by exploring several application problems that employ a combination of mesh refinement, mesh motion, and temporal refinement.

I will percent analytical and numerical results of the PFASST algoritm applied to something other than the heat equation.

Spectral methods solve elliptic PDEs numerically with errors bounded by an exponentially decaying function of the number of modes when the solution is analytic. For time dependent problems, almost all focus has been on low-order finite difference schemes for the time derivative and spectral schemes for spatial derivatives. Spectral methods which converge spectrally in both space and time have appeared recently. This paper shows the exponential convergence of the heat equation for a Legendre spectral collocation method. A condition number estimate of the method is also given.

Spectral deferred corrections (SDC) are an easy way to construct higher-order time integration schemes from simple base integrators, e.g. backward Euler or velocity-Verlet. In addition, their iterative nature makes them an attractive subject for investigating parallelization in time. The most prominent example is the "parallel full approximation scheme in space in time" (PFASST), which couples SDC and multigrid techniques in a clever yet slightly intricate way. While PFASST parallelizes the computation of multiple steps in time, we will focus on the direct, single-step parallelization of SDC iterations themselves. In this talk, we will explore different approaches and show their strengths (if any) as well as their (often severe) limitations. Results of this investigation will turn out to be rather diverse, so this talk is also meant to challenge the community to find better and more robust ways for parallelizing SDC across the method.

The parareal in time algorithm, combines a serial coarse solver used on the full time simulation with loops of fine solver, implemented in parallel, over shorter time simulations. This algorithm uses a new direction for task decomposition and parallelism for time dependent (partial) differential equation. One way to improve the parallel efficiency of the parareal in time algorithm, is to degrade the fist swap of the fine solver and improve its accuracy over the parallel iterations. This degradation and improvement over the loops has to be well balanced in order not to pollute the convergence efficiency of the method. This presentation will propose some mathematical approach for rigorously sustain this approach.

Since the introduction of the idea of adding parallelism to time integration in the seminal work of Nievergelt in 1964, various approaches for parallel-in-time integration have been explored. One of these approaches is applying multigrid to the time dimension, resulting in the multigrid-reduction-in-time (MGRIT) algorithm by Falgout et al. The corresponding open source code XBraid is flexible, allowing for a variety of time stepping, relaxation, and temporal coarsening options. It was recently shown in Bolten et al. that the time-parallel expansion of spectral deferred corrections methods, the Parallel Full Approximation Scheme in Space and Time (PFASST) algorithm by Emmett and Minion, can be described as a multigrid-in-time method under certain assumptions. In this talk, we discuss approaches for incorporating aspects of PFASST and multigrid-based parallel-in-time methods into XBraid.

Communication bandwidth between nodes of parallel machines grows slower than the nodes' computational performance and increasingly becomes a bottleneck of PDE solvers. We present lossy compression techniques based on multilevel prediction in hierarchical grids and entropy coding as a means to reduce the amount of data to be communicated between nodes in SDC-based time-parallel simulations. Compression rate and impact on convergence depending on the quantization error are investigated both theoretically and at some numerical examples realized in the finite element code Kaskade 7, which leads to an adaptive choice of the quantization threshold.

For the model problem of the heat equation we formulate and describe space-time finite and boundary element methods. Here the discretization is done with respect to rather general discretizations of the space-time cylinder and its boundary. This allows the use of adaptivity both in space and time, but requires the solution of the global system at once. Hence, preconditioning and parallelization are mandatory. Numerical results are given, including first results for the wave equation.

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.