David C. Seal

Andrew J. Christieb, Sigal Gottlieb, Zachary J. Grant et al.

High order strong stability preserving (SSP) time discretizations are advantageous for use with spatial discretizations with nonlinear stability properties for the solution of hyperbolic PDEs. The search for high order strong stability time-stepping methods with large allowable strong stability time-step has been an active area of research over the last two decades. Recently, multiderivative time-stepping methods have been implemented with hyperbolic PDEs. In this work we describe sufficient conditions for a two-derivative multistage method to be SSP, and find some optimal SSP multistage two-derivative methods. While explicit SSP Runge--Kutta methods exist only up to fourth order, we show that this order barrier is broken for explicit multi-stage two-derivative methods by designing a three stage fifth order SSP method. These methods are tested on simple scalar PDEs to verify the order of convergence, and demonstrate the need for the SSP condition and the sharpness of the SSP time-step in many cases.

Andrew J. Christlieb, Xiao Feng, David C. Seal et al.

We propose a high-order finite difference weighted ENO (WENO) method for the ideal magnetohydrodynamics (MHD) equations. The proposed method is single-stage, single-step, maintains a discrete divergence-free condition on the magnetic field, and has the capacity to preserve the positivity of the density and pressure. To accomplish this, we use a Taylor discretization of the Picard integral formulation (PIF) of the finite difference WENO method proposed in [SINUM, 53 (2015), pp. 1833--1856], where the focus is on a high-order discretization of the fluxes (as opposed to the conserved variables). We use the version where fluxes are expanded to third-order accuracy in time, and for the fluid variables space is discretized using the classical fifth-order finite difference WENO discretization. We use constrained transport in order to obtain divergence-free magnetic fields, which means that we simultaneously evolve the magnetohydrodynamic and magnetic potential equations, and set the magnetic field to be the (discrete) curl of the magnetic potential after each time step. In this work, we compute these curls to fourth-order accuracy. In order to retain a single-stage, single-step method, we develop a novel Lax-Wendroff discretization for the evolution of the magnetic potential, where we start with technology used for Hamilton-Jacobi equations in order to construct a non-oscillatory magnetic field. Positivity preservation is realized by introducing a parameterized flux limiter that considers a linear combination of high and low-order numerical fluxes. This positivity limiter lacks energy conservation. However, this limiter can be dropped for problems where the pressure does not become negative. We present two and three dimensional numerical results for several standard test problems. These results assert the robustness and verify the high-order of accuracy of the proposed scheme.

Matthew F. Causley, Hana Cho, Andrew J. Christlieb et al.

We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy in space and time. The solver is built by first constructing a single-dimensional heat equation solver that uses fast O(N) convolution. This fundamental solver has arbitrary order of accuracy in space, and is based on the use of the Green's function to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and convergence, and permit us to construct schemes that have provable stiff decay. The multi-dimensional solver is built by repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and make use of Hermite-Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear and nonlinear equations including heat, Allen-Cahn, and the Fitzhugh-Nagumo system of equations in one and two dimensions.

David C. Seal, Qi Tang, Zhengfu Xu et al.

In this work we construct a high-order, single-stage, single-step positivity-preserving method for the compressible Euler equations. Space is discretized with the finite difference weighted essentially non-oscillatory (WENO) method. Time is discretized through a Lax-Wendroff procedure that is constructed from the Picard integral formulation (PIF) of the partial differential equation. The method can be viewed as a modified flux approach, where a linear combination of a low- and high-order flux defines the numerical flux used for a single-step update. The coefficients of the linear combination are constructed by solving a simple optimization problem at each time step. The high-order flux itself is constructed through the use of Taylor series and the Cauchy-Kowalewski procedure that incorporates higher-order terms. Numerical results in one- and two-dimensions are presented.

Scott A. Moe, James A. Rossmanith, David C. Seal

This work introduces a single-stage, single-step method for the compressible Euler equations that is provably positivity-preserving and can be applied on both Cartesian and unstructured meshes. This method is the first case of a single-stage, single-step method that is simultaneously high-order, positivity-preserving, and operates on unstructured meshes. Time-stepping is accomplished via the Lax-Wendroff approach, which is also sometimes called the Cauchy-Kovalevskaya procedure, where temporal derivatives in a Taylor series in time are exchanged for spatial derivatives. The Lax-Wendroff discontinuous Galerkin (LxW-DG) method developed in this work is formulated so that it looks like a forward Euler update but with a high-order time-extrapolated flux. In particular, the numerical flux used in this work is a linear combination of a low-order positivity-preserving contribution and a high-order component that can be damped to enforce positivity of the cell averages for the density and pressure for each time step. In addition to this flux limiter, a moment limiter is applied that forces positivity of the solution at finitely many quadrature points within each cell. The combination of the flux limiter and the moment limiter guarantees positivity of the cell averages from one time-step to the next. Finally, a simple shock capturing limiter that uses the same basic technology as the moment limiter is introduced in order to obtain non-oscillatory results. The resulting scheme can be extended to arbitrary order without increasing the size of the effective stencil. We present numerical results in one and two space dimensions that demonstrate the robustness of the proposed scheme.

Mathew F. Causley, David C. Seal

In this work we analyze the convergence properties of the Spectral Deferred Correction (SDC) method originally proposed by Dutt et al. [BIT, 40 (2000), pp. 241--266]. The framework for this high-order ordinary differential equation (ODE) solver is typically described wherein a low-order approximation (such as forward or backward Euler) is lifted to higher order accuracy by applying the same low-order method to an error equation and then adding in the resulting defect to correct the solution. Our focus is not on solving the error equation to increase the order of accuracy, but on rewriting the solver as an iterative Picard integral equation solver. In doing so, our chief finding is that it is not the low-order solver that picks up the order of accuracy with each correction, but it is the underlying quadrature rule of the right hand side function that is solely responsible for picking up additional orders of accuracy. Our proofs point to a total of three sources of errors that SDC methods carry: the error at the current time point, the error from the previous iterate, and the numerical integration error that comes from the total number of quadrature nodes used for integration. The second of these two sources of errors is what separates SDC methods from Picard integral equation methods; our findings indicate that as long as difference between the current and previous iterate always gets multiplied by at least a constant multiple of the time step size, then high-order accuracy can be found even if the underlying "solver" is inconsistent the underlying ODE. From this vantage, we solidify the prospects of extending spectral deferred correction methods to a larger class of solvers to which we present some examples.

Jochen Schütz, David C. Seal, Alexander Jaust

In this work, we construct novel discretizations for the unsteady convection-diffusion equation. Our discretization relies on multiderivative time integrators together with a novel discretization that reduces the total number of unknowns for the solver. These type of temporal discretizations come from an umbrella class of methods that include Lax-Wendroff (Taylor) as well as Runge-Kutta methods as special cases. We include two-point collocation methods with multiple time derivatives as well as a sixth-order fully implicit collocation method that only requires a total of three stages. Numerical results for a number of sample linear problems indicate the expected order of accuracy and indicate we can take arbitrarily large time steps.

Alexander Jaust, Jochen Schütz, David C. Seal

In this paper we apply implicit two-derivative multistage time integrators to viscous conservation laws in one and two dimensions. The one dimensional solver discretizes space with the classical discontinuous Galerkin (DG) method, and the two dimensional solver uses a hybridized discontinuous Galerkin (HDG) spatial discretization for efficiency. We propose methods that permit us to construct implicit solvers using each of these spatial discretizations, wherein a chief difficulty is how to handle the higher derivatives in time. The end result is that the multiderivative time integrator allows us to obtain high-order accuracy in time while keeping the number of implicit stages at a minimum. We show numerical results validating and comparing methods.

David C. Seal, Yaman Güçlü, Andrew J. Christlieb

High-order temporal discretizations for hyperbolic conservation laws have historically been formulated as either a method of lines (MOL) or a Lax-Wendroff method. In the MOL viewpoint, the partial differential equation is treated as a large system of ordinary differential equations (ODEs), where an ODE tailored time-integrator is applied. In contrast, Lax-Wendroff discretizations immediately convert Taylor series in time to discrete spatial derivatives. In this work, we propose the Picard integral formulation (PIF), which is based on the method of modified fluxes, and is used to derive new Taylor and Runge-Kutta (RK) methods. In particular, we construct a new class of conservative finite difference methods by applying WENO reconstructions to the so-called "time-averaged" fluxes. Our schemes are automatically conservative under any modification of the fluxes, which is attributed to the fact that classical WENO reconstructions conserve mass when coupled with forward Euler time steps. The proposed Lax-Wendroff discretization is constructed by taking Taylor series of the flux function as opposed to Taylor series of the conserved variables. The RK discretization differs from classical MOL formulations because we apply WENO reconstructions to time-averaged fluxes rather than taking linear combinations of spatial derivatives of the flux. In both cases, we only need one projection onto the characteristic variables per time step. The PIF is generic, and lends itself to a multitude of options for further investigation. At present, we present two canonical examples: one based on Taylor, and the other based on the classical RK method. Stability analyses are presented for each method. The proposed schemes are applied to hyperbolic conservation laws in one- and two-dimensions and the results are in good agreement with current state of the art methods.

Scott A. Moe, James A. Rossmanith, David C. Seal

The discontinuous Galerkin (DG) finite element method when applied to hyperbolic conservation laws requires the use of shock-capturing limiters in order to suppress unphysical oscillations near large solution gradients. In this work we develop a novel shock-capturing limiter that combines key ideas from the limiter of Barth and Jespersen [AIAA-89-0366 (1989)] and the maximum principle preserving (MPP) framework of Zhang and Shu [Proc. R. Soc. A, 467 (2011), pp. 2752--2776]. The limiting strategy is based on traversing the mesh element-by-element in order to (1) find local upper and lower bounds on user-defined variables by sampling these variables on neighboring elements, and (2) to then enforce these local bounds by minimally damping the high-order corrections. The main advantages of this limiting strategy is that it is simple to implement, effective at shock capturing, and retains high-order accuracy of the solution in smooth regimes. The resulting numerical scheme is applied to several standard numerical tests in both one and two-dimensions and on both Cartesian and unstructured grids. These tests are used as benchmarks to verify and assess the accuracy and robustness of the method.