Computer simulations of evolutionary multiscale multiphysics partial differential equations are important in many areas of science and engineering. Algorithms for time integration of these systems face important challenges. Multiscale problems have components evolving at different rates. No single time step can solve all components efficiently (e.g., when an explicit discretization is used, and the spatial discretization uses both fine and coarse mesh patches). Multiphysics problems are driven by multiple simultaneous processes with different dynamic characteristics. No single time discretization method is best suited to solve all processes (e.g., when some are stiff and others non-stiff).
In order to address these challenges, multimethods have been proposed. Multimethods are time integration approaches that use different solution strategies for different subsystems have been developed. For example, different processes are discretized with different numerical schemes, and different components of the system are solved with different time steps. We discuss several general aspects of multimethods for the integration for multiphysics systems, as well as new developments in the field.

The dispersive shallow water equations (DSWE) are fluid models, applicable to coastal regions that include additional physics (such as dispersion) to the well-known shallow water equations. The DSWEs present several challenges for efficient time-stepping including mixed space and time derivatives, nonlinearities and higher order spatial derivatives. We present semi-implicit high-order time stepping strategies that avoid a fully implicit treatment of the nonlinear terms and simplify the treatment of the mixed space-time derivatives. The approach is based on extensions of a recent unconditional stability theory, where due to the structure of the equations, zero stability plays the role of unconditional stability.

This talk introduces new partitioned time integration methods aimed at practical PDE problems with different dynamical characteristics. The proposed methods allot integrators with different properties to different PDE components to reduce the time-to-solution for a desired accuracy. Classical strategies such as IMEX and multirate will be combined to achieve this purpose under such constraints as conservation of linear invariants and entropy stability. A strong emphasis is placed on what practitioners identify as essential components to be resolved, influencing how one interprets accuracy.

In recent years computational simulations have rapidly evolved in complexity (high order discretizations, spatial adaptivity, additional physical processes), placing ever larger strains on the time integration methods on which they rely. High spatial order necessitates comparably high order time integration. Spatial adaptivity and multiphysics processes give rise to subsets of the solution that evolve at differing time scales, or to simulations that combine nonstiff but highly nonlinear processes with others that may be highly stiff but that are frequently linear. In this talk, I will discuss recent work on advanced time integration methods (both mathematical derivation and implementation in high-quality open-source software), that allow extreme flexibility in the use of different techniques to distinct physical processes, while still allowing high orders of accuracy. I will primarily focus on the newly-developed IMEX-MRI-GARK methods (Chinomona and R., 2021), and their implementation in the ARKODE library within SUNDIALS, but I will additionally point out other recent related work.

The ASPECT code (short for "Advanced Solver for Problems in Earth ConvecTion") solves problems relating to the motion of material, energy, and composition in both the Earth crust and the mantle. This is described by a coupled system that involves instantaneous Stokes flow, the advection and diffusion of a temperature field, and the advection and reaction of chemical compositions or other quantities. It can also simulateously solve for the deformation of the Earth surface, as well as advect particles with the flow field. The resulting time stepping scheme is diverse, but for the most part only uses home-grown time steppers that, for historical reasons, are not well integrated.
I will discuss what we do and why, as well as what we do not do and why not. I will also discuss what I think better solutions could be.

Time filtering has been used to enhance the order of accuracy of given methods. This is particularly useful in the context of legacy codes, in which the time-stepping module is given and difficult to change. However, modifying the inputs and outputs is simple and allows for higher order. In this talk, we show how time filtering approaches can be seen as equivalent to generating a general linear methods. We use this GLM approach to develop an optimization routine that allows us to find new time-filtering methods with high order and efficient linear stability properties. We will present our new methods and show their performance when tested on sample problems.

Authors: Carol S. Woodward(1), Daniel R. Reynolds(2), David J. Gardner(1), and Cody J. Balos(1)
Affiliations: (1) Lawrence Livermore National Laboratory, (2) Southern Methodist University
As the scientific computing community leverages unprecedented high-performance computing resources to work toward predictive science, software complexity is increasing. Applications increasingly require newly developed numerical methods implemented for a variety of programming environments as well as the combined use of independent software packages, which have diverse sponsors, priorities, and processes for development and release. Numerical software libraries provide a large and growing resource for high-quality, reusable software components on which applications can be rapidly constructed, with improved robustness, portability, and sustainability.
SUNDIALS is a suite of robust and scalable integrators and solvers for systems of ordinary differential equations, differential-algebraic equations, and nonlinear equations designed for use on computing systems ranging from desktop machines to supercomputers. 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 SUNDIALS history and capabilities. We will also discuss the design principles adopted by the SUNDIALS development team and discuss how they are manifested in package flexibility and user interfaces. We will give examples of how this flexibility provides an insertion path for new methods developed in the community to get into applications. This presentation will also introduce the xSDK, or Extreme-scale Scientific Software Development Kit, where community-defined policies are increasing the quality and interoperability across numerical libraries as needed by the DOE Exascale Computing Project. This talk will present the community policies developed through the xSDK project, their expected benefits, and lessons learned from adoption of these policies into existing packages. The presentation will conclude with examples highlighting applications of SUNDIALS in some DOE simulations and give a preview of new developments.
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- 830479.

Over the past decades, exponential integration emerged as a numerical technique that can offer significant computational advantages for integrating large-scale stiff systems of differential equations. In this talk, we will provide an overview of the ways exponential approach can be used to construct new methods and discuss theoretical and practical aspects of design, analysis and implementation of different classes of exponential integrators. We will illustrate performance gains these schemes provide using test problems and several applications in fluid mechanics, plasma physics and computer graphics.

In the 2007 Sidney Fernbach Award lecture, I argued for going fully implicit in time-dependent simulations in which the stability restriction for explicit methods is more stringent than is required for accuracy in resolving the phenomena of interest. These include multirate problems with good scale separation, coupled multiphysics problems, and “outer loop” problems in design, control, inversion, and assimilation. The basis for the argument was that scalable and effective algorithms are known and available in open-source software that harvest much legacy code in “call-back” form and thus offer a reasonable learning curve. While none of this has disappeared, much new has appeared that forces consideration of additional algorithmic and software infrastructure in bridging the widening gap between ambitious applications and austere architectures. We comment on the evolution of this infrastructure and attempt to integrate approaches described in the workshop under a bigger tent.

The efficient use of modern supercomputers has become one of the key challenges in computational science. For the numerical solution of time-dependent processes, time-parallel methods have opened new ways to overcome both strong and weak scaling limits. If higher order accuracy in time is needed or at least feasible, parallelization techniques based on collocation methods can provide temporal parallelism within a single time-step as well as across multiple time-steps. In this talk, we give an overview of three different approaches to introduce parallel-in-time integration for collocation methods: (1) parallel preconditioners for spectral deferred corrections, (2) the parallel full approximation scheme in space and time, and (3) diagonalization-based preconditioning for multi-step collocation problems. Those approaches can even be combined to obtain multi-time-parallel integrators. We shed light on the pros and cons of the different variants, their implementation on HPC systems, as well as existing and potential applications. We also discuss current roadblocks and further research directions.

We discuss time stepping challenges in large scale PDE solvers built within an adaptive meshing framework. Adaptive meshing is a software strategy for dynamically managing mesh resolution so that spatial and temporal solution features of interest (typically defined by the user) are computed at the highest levels of resolution. By dynamically managing computational resources, solver performance can be dramatically increased. Most commonly, only the spatial resolution is dynamically managed. In this case, a standard "method of lines" approach can be used to advance the solution, with possible time step restrictions based on the highest spatial resolution. Less common is the use of "local time stepping", a strategy in which the time step size is locally adjusted to maintain, for example, a constant CFL number across the domain. Local time stepping works well for purely hyperbolic problems, and in fact is implemented in many Cartesian grid codes. What is not entirely clearly is how to manage local time stepping with multi-stage time stepping schemes or operator split approaches for systems of PDEs, especially if the system involves elliptic terms.
I will discuss these issues in the context of ForestClaw software (www.forestclaw.org), a Cartesian grid based software platform for solving PDEs on a hierarchy of adaptively refined Cartesian meshes. We will show progress on several projects, including a DARPA funded project involving remote sensing in the atmosphere. This project will couple several complex codes through the ForestClaw platform. The coupling will occur through interpolation between different ForestClaw mesh hierarchies, as well as through operator splitting involving elliptic, parabolic and hyperbolic terms. Time stepping must be commensurate with the accuracy of each individual code, while balancing communication costs.

Upcoming exascale architectures require rethinking of the numerical algorithms used in large-scale PDE-based applications. These architectures favor algorithms, such as high-order finite elements, that expose fine-grain parallelism and maximize the ratio of floating point operations to energy intensive data movement. In this talk we present an overview of MFEM [1], a scalable library for high-order finite element discretization of PDEs on general unstructured grids. We also report on recent work in the Center for Efficient Exascale Discretizations [2], a co-design center in the US Exascale Computing Project focused on next-generation discretization software and algorithms. We will describe recent research on performance optimizations for GPU architectures, scalable unstructured adaptive mesh refinement and matrix-free preconditioning, with special attention to the needs for time integration methods in high-order applications.

[1] MFEM: Modular finite element library, http://mfem.org.
[2] Center for Efficient Exascale Discretizations, http://ceed.exascaleproject.org.

In this talk, we overview a class of high order conservative semi-Lagrangian schemes for BGK model. The schemes are constructed by coupling the conservative non-oscillatory reconstruction technique with a conservative treatment of the collision term, obtained by either a discrete Maxwellian or by an L2-minimization technique. Due to the semi-Lagrangian nature, the time step is not restricted by a CFL-type condition, while the implicit treatment of the relaxation term based on Runge-Kutta or BDF time discretization enables us to avoid the stiffness problem coming from a small Knudsen number. Because of L-stability and exact conservation, the resulting scheme is asymptotic preserving for the underlying fluid dynamic limit. Several test cases confirm the accuracy and robustness of the methods, and the AP property of the schemes. The method has been extended to the treatment of inert gas mixtures, and applied to compare different models in various regimes.
In general, such approaches use fixed velocity grids, and one must secure a sufficient number of grid points in phase space to resolve the structure of the distribution function. When dealing with high Mach number problems, where large variation of mean velocity and temperature are present in the domain under consideration, the computational cost and memory allocation requirements become prohibitively large. Local velocity grid methods have been developed to overcome such difficulty in the context of Eulerian based schemes. In this talk, we introduce a velocity adaption technique for the semi-Lagrangian scheme applied to the BGK model. The velocity grids will be set locally in time and space. We apply a weighted minimization approach to impose global conservation, generalizing the L2-minimization technique introduced in. The efficiency of the proposed scheme is illustrated numerically.
An additional application of conservative SL schemes concerns the numerical simulation of Vlasov-type equations. Here conservation of the scheme will provide some advantage over standard non conservative schemes for long time computation. The research is conducted with the following collaborators: S. Boscarino, S. Y. Cho, M. Groppi, S.-B. Yun, JM Qiu, and T. Xiong.

Order reduction, i.e., the convergence of the solution at a lower rate than the formal order of the chosen time-stepping scheme, is a fundamental problem in stiff ODEs, and particularly in PDE IBVPs. Runge-Kutta schemes with high stage order provide a remedy, but unfortunately high stage order is incompatible with DIRK schemes. We first highlight the spatial manifestations of order reduction in PDE IBVPs. Then we introduce the concept of weak stage order, and (a) demonstrate how it overcomes order reduction in important linear PDE problems; and (b) how high-order DIRK schemes can be constructed that are devoid of order reduction.

We propose an Eulerian-Lagrangian (EL) Runge-Kutta (RK) discontinuous Galerkin (DG) method. The method is designed as a generalization of the semi-Lagrangian (SL) DG method, which is formulated based on an adjoint problem and tracing upstream cells by tracking characteristics curves highly accurately. Depending on the velocity field, the shape of upstream cells could be of arbitrary shape, for which a more sophisticated approximation is required to get high order approximation. For example, in the SLDG algorithm, quadratic-curved (QC) quadrilaterals were proposed to approximate upstream cells in order to obtain third order spatial accuracy in a swirling deformation example. In this work, we propose a more general formulation, named the ELDG method, for which the scheme is based on a modified adjoint problem for which upstream cells are always quadrilaterals. This leads to a new formulation of ELDG method, which avoids the need to use QC quadrilaterals to better approximate upstream cells in the original SLDG algorithm. The newly proposed ELDG method can be viewed as a new general framework, in which both the classical RK DG formulation and the SL DG formulation can fit in. Numerical results on linear advection problems, as well as the nonlinear Vlasov dynamics using the exponential RK framework, will be presented to demonstrate the effectiveness of the proposed approach.

We present a framework for compatible spatial discretizations and time integrators that incorporates properties of the combined space-time integrals of a given PDE. Demonstrations for model equations will be provided, ranging from hyperbolic conservation laws, such as Lagrangian advection and Maxwell’s equations, and parabolic systems like heat transfer with phase change. The approach is based on a “cut-cell” discretization that tracks interfaces and jump conditions of the PDE along an approximate space-time boundary. It can also be extended to include moving meshes, and even mesh discontinuities in space-time, which otherwise would introduce large errors and overly-restrictive time steps for stability with typical approaches. We'll demonstrate that in comparison with a generic method-of-lines discretization, better accuracy and stability properties can be achieved when considering design of space- and time-discretizations together.

In this work we present a class of high order unconditionally strong stability preserving (SSP) implicit multi-derivative Runge--Kutta schemes, and SSP implicit-explicit (IMEX) multi-derivative Runge--Kutta schemes where the time-step restriction is independent of the stiff term. The unconditional SSP property for a method of order p>2 is unique among SSP methods, and depends on a backward-in-time assumption on the derivative of the operator. We show that this backward derivative condition is satisfied in many relevant cases where SSP IMEX schemes are desired. We devise unconditionally SSP implicit Runge--Kutta schemes of order up to p=4, and IMEX Runge--Kutta schemes of order up to p=3. For the multi-derivative IMEX schemes, we also derive and present the order conditions, which have not appeared previously. The unconditional SSP condition ensures that these methods are positivity preserving, and we present sufficient conditions under which such methods are also asymptotic preserving when applied to a range of problems, including a hyperbolic relaxation system, the Broadwell model, and the Bhatnagar-Gross-Krook (BGK) kinetic equation. We present numerical results to support the theoretical results, on a variety of problems. This is joint work with S. Gottlieb, Z. Grant, and R. Shu.

In this talk, I will review the evolution of the summation-by-parts (SBP) framework. Starting from linear PDEs, I will discuss how this framework has matured from its finite-difference origins into an abstract matrix analysis framework that is nearly discretization agnostic for the discretization of spacial derivatives. I will next summarize work by several groups on the link between implicit SBP operators for temporal derivatives and Runge-Kutta (RK) methods and pose a number of questions for discussion. The compelling features of the SBP framework are that it enables the analysis and modification of the actual algorithms implemented in practice (e.g., it accounts for variational crimes such as inexact integration) and leads to the construction of schemes with provable properties (e.g., stability and conservation). I will then move to nonlinear conservation laws where at the continuous level stability can be proven via entropy-stability analysis and demonstrate how these same ideas (and stability proofs) can be constructed leveraging the SBP framework and Tadmor’s two-point flux functions for both space and time. I will then cover alternatives to using implicit SBP operators for temporal derivatives in the form of RK relaxation schemes and (time permitting) a brief review of new work by Yamaleev and Upperman on constructing positivity preserving schemes. I will finish with the nonlinear analogue of the questions I presented in the linear section.

Rothe’s method, transvers method of lines and the Method of Lines Transpose are all approaches for solving time dependent PDEs that approach the problem by turning the PDE into a boundary value problem and then addressing the resulting BVP with an efficient spatial discretization. Traditionally, this approach was coupled with an implicit formulation of a time marching method as the very first step in setting up the numerical integrator. In the Method of Lines Transpose approach, the resulting BVP has been addressed using kernel based methods and efficiently evaluated using fast kernel tricks. In this talk, we discuss a refactoring of these kernel based methods into a form where differential operators of any PDE can be express through Successive Convolution independent of the time integration strategy. The method leads to a formulation that makes explicit time stepping methods provably unconditionally stable for linear problems and behaves unconditionally stable for non-linear problems. The method is O(N) and relies on WENO integration to address problems with discontinuities. Boundary conditions pose an interesting challenge for the method, and will be discussed during this talk. We demonstrate the method by a applying it to the Hamilton Jacobi and Degenerate Advection Diffusion equations in 1D and 2D.
This is joint work with my student, Mr. Bill Sands, and three former post docs, Drs. Hyoseon Yang, Yan Jiang and Wei Gou.

Many mathematical models are equipped with an energy that is conserved or an entropy that is known to change monotonically in time. Integrators that preserve these properties discretely are usually expensive, with the best-known examples being fully-implicit Runge-Kutta methods. I will present a modification that can be applied to any integrator in order to preserve such a structural property. The resulting method can be fully explicit, or (depending on the functional) may require the solution of a scalar algebraic equation at each step. I will present examples to show the effectiveness of these “relaxation” methods, and their advantages over fully implicit methods or orthogonal projection. Examples will include applications to compressible fluid dynamics, dispersive nonlinear waves, and Hamiltonian systems.