Strong Stability Preserving Integrating Factor Runge-Kutta Methods Sigal Gottlieb Joint work with: Leah Isherwood and Zachary J. Grant Advances in PDEs: Theory, Computation and Application to CFD ICERM August 20-24, 2018 Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 1 / 38
Overview In memory of Saul (Shalom) Abarbanel 1 Strong Stability Preserving (SSP) 2 Strong Stability Preserving Runge-Kutta (SSPRK) methods Exponential Integrators Motivation 3 Integrating Factor (IF) methods Integrating Factor Runge-Kutta (IFRK) methods Strong Stability Preserving Integrating Factor Runge-Kutta (SSPIFRK) methods Numerical Results 4 Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 2 / 38
Transformations: Saul/Shalom Communication: Joan and Bob Optimality: Postdoc positions Relativity: Uncle Shalom to me and my brothers Uniqueness: Many Friday night dinners over 40+ years Positivity preservation: ICASE names, insult-work sessions Energy maximization Travel and NSF panels Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 3 / 38
Strong Stability Preserving (SSP) Motivation Consider the hyperbolic partial differential equation U t + f ( U ) x = 0 . Method of lines approach: we discretize the problem in space, to obtain some ODE of the form u t = F ( u ) and we evolve this ODE in time using standard time-stepping methods such as Runge–Kutta methods. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 4 / 38
Strong Stability Preserving (SSP) Motivation Given a linear differential equation and consistent linear numerical method, the Lax-Richtmyer equivalence theorem tells us that linear stability is necessary and sufficient for convergence. consistency + stability ⇐ ⇒ convergence For nonlinear PDEs, if a numerical method is consistent and its linearization is L 2 stable and adequately dissipative, then for sufficiently smooth problems the nonlinear approximation is convergent (Strang 1964). ⇒ For smooth solutions, we look at L 2 linear stability (plus some dissipativity) to prove convergence. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 5 / 38
Strong Stability Preserving (SSP) Motivation Hyperbolic conservation laws tend to have solutions which start with or develop discontinuities or steep gradients over time. 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.6 −0.8 −0.8 −1 −1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 If the solution is discontinuous , then we can get oscillations, non-physical negative values, etc. For a nonlinear problem with discontinuous solutions, linear L 2 stability analysis is is not sufficient for convergence. We build spatial discretizations which satisfy some nonlinear, non-inner product stability properties. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 6 / 38
Strong Stability Preserving (SSP) Consider ODE system u t = F ( u ) , where spatial discretization F ( u ) is carefully chosen, (e.g. TVD, TVB, ENO, WENO, positivity or maximum principle preserving) so that the solution from the forward Euler method u n +1 = u n + ∆ tF ( u n ) , satisfies the monotonicity (or strong stability) requirement || u n +1 || ≤ || u n || , in some norm, semi-norm, or convex functional || · || , for a suitably restricted timestep ∆ t ≤ ∆ t FE . Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 7 / 38
Strong Stability Preserving (SSP) Motivation We now have a spatial discretization that is stable when coupled with forward Euler . But in practice, we need to use higher order methods. There has been much work designing spatial discretizations which satisfy certain nonlinear , non-inner product stability properties when coupled with Forward Euler. Forward Euler is only first order accurate error= O (∆ t ) Linear stability region fails to capture the imaginary axis These issues can be handled by using a higher order time integrator. How can these time integrators also preserve the monotonicity property guaranteed by the Forward Euler time step? Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 8 / 38
Strong Stability Preserving (SSP) Given an operator F ( u ) with forward Euler monotonicity condition � u n +1 � = � u n + ∆ tF ( u n ) � ≤ � u n � under time-step restriction ∆ t ≤ ∆ t FE . Creating the higher order time integrator is done by Decomposing a higher order time-stepping method into convex combinations of forward Euler so that any convex-functional monotonicity property � u n +1 � ≤ � u n � will be preserved under a time step restriction ∆ t ≤ C ∆ t FE . Any higher order method that can be decomposed in this way is called strong stability preserving with SSP coefficient C . This convex combination condition is both necessary and sufficient. Decoupling the analysis C is a property only of the time-integrator. ∆ t FE is a property of the spatial discretization. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 9 / 38
Example: Shu-Osher third order method For example, Shu-Osher’s third order SSP method: eSSPRK(3,3): u n + ∆ tF ( u n ) u (1) = 3 4 u n + 1 � u (1) + ∆ tF ( u (1) ) � u (2) = 4 1 3 u n + 2 � u (2) + ∆ tF ( u (2) ) � u n +1 = 3 This three stage method has an SSP coefficient C =1. i.e ∆ t = ∆ t FE The notation eSSPRK(s,p) denotes an explicit SSP Runge–Kutta method with s stages and of order p . Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 10 / 38
SSP Runge-Kutta methods For example, an s -stage explicit Runge–Kutta method can be written as: u (0) u n , = i − 1 α i , j u ( j ) + ∆ t β i , j F ( u ( j ) ) � � � u ( i ) = , i = 1 , ..., s j =0 u n +1 u ( s ) . = If all the coefficients α i , j and β i , j are non-negative, and a given α i , j is zero only if its corresponding β i , j is zero,then each stage can be rearranged into a convex combination of forward Euler steps � � i − 1 i − 1 � � � � � u ( j ) + ∆ t β i , j � α i , j u ( j ) + ∆ t β i , j F ( u ( j ) ) � � � � u ( i ) � = � � � F ( u ( j ) ) � ≤ � u n � � ≤ α i , j � � � � α i , j � � j =0 j =0 � � provided that the time-step satisfies α i , j ∆ t ≤ min ∆ t FE . β i , j i , j Note that for consistency � i − 1 j =0 α i , j = 1. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 11 / 38
SSP Runge–Kutta Methods: bounds and barriers Explicit SSP Runge–Kutta methods Order barrier p ≤ 4 SSP coefficent C ≤ s Implicit SSP Runge–Kutta methods Order barrier p ≤ 6 SSP coefficient C ≤ 2 s for order p ≥ 2 Implicit-explicit (IMEX) SSP Runge–Kutta methods Order barrier p ≤ 4 SSP coefficient between C ≤ s and C ≤ 2 s Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 12 / 38
Exponential Integrators Motivation We consider a problem of the form u t = Lu + N ( u ) As before, the two components satisfy � u n + ∆ tN ( u n ) � ≤ � u n � for ∆ t ≤ ∆ t FE while taking a forward Euler step using the linear component Lu results in the strong stability condition � u n + ∆ tLu n � ≤ � u n � ∆ t ≤ ˜ for ∆ t FE where ˜ ∆ t FE << ∆ t FE . Our focus is how to preserve the nonlinear, non-inner product monotonicity properties, while avoiding the severe time-step restriction coming from the linear component. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 13 / 38
Motivating example Consider � � 1 � 1 , if 0 ≤ x ≤ 1 / 2 2 u 2 u t + 10 u x + = 0 u (0 , x ) = 0 , if x > 1 / 2 x using first order upwind finite difference. The relevant TVD time-step restrictions are: Explicit SSPRK(3,3): ∆ t ≤ 0 . 09∆ x Implicit SSPRK(2,3): ∆ t ≤ 0 . 248∆ x IMEX SSP(3,3, 1 10 ): ∆ t ≤ 0 . 149∆ x The time-step restriction is driven by the linear component – even when handled implicitly! When linear L 2 stability is the concern, implicit or IMEX methods completely resolve the stiffness coming from the linear case. However, when we require the SSP condition to hold, implicit or IMEX do not significantly alleviate the time-step restriction. Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 14 / 38
Exponential Integrators background Emerged in the 1960s as an alternative to implicit methods for stiff problems Rather than solving large linear systems these methods require the exponential of a matrix (typically linear operator or Jacobian) Using direct methods to evaluate the matrix exponential such as Taylor or Pad´ e approximations Direct approximations made exponential integrator methods too costly for large scale problems Focus was on small scale problems Re-emerged in 1980s Using iterative techniques to compute the matrix exponential Enables exponential integrator methods to be used on large scale problems Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 15 / 38
Exponential Integrators Starting from the equation u t = Lu + N ( u ) we solve the linear component by an integrating factor e − Lt u t − e − Lt Lu = e − Lt N ( u ) � � e − Lt u t = e − Lt N ( u ) . The following classes of exponential integrators Integrating factor methods (IF) 1 Exponential time differencing methods (ETD) differ in how they approximate the nonlinear component. 1 Also known as Lawson methods Sigal Gottlieb (UMassD) SSPIFRK ICERM 2018 16 / 38
Recommend
More recommend
