Designing Time Filters Using a General Linear Method Approach Sigal - PowerPoint PPT Presentation

Designing Time Filters Using a General Linear Method Approach Sigal Gottlieb (UMassD) and Zachary J. Grant (ORNL) Joint work with: Victor P. DeCaria & William J. Layton (UPitt) International Conference on Advances in Applied Mathematics at Tel Aviv University December 18-20, 2018 Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 1 / 26

Overview In memory of Saul (Shalom) Abarbanel 1 Background: Time Filters 2 Time Filters as General Linear Methods 3 Methods found using our GLM approach and optimization code 4 Conclusions 5 Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 2 / 26

Memories of Shalom Uncle Shalom to me and my brothers Many Friday night dinners over 40+ years Travel and NSF panels Advice Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 3 / 26

Motivation Consider Layton’s filtered implicit Euler scheme: u n + ∆ tF ( y (1) ) y (1) = y (1) − 1 � u n − 1 − 2 u n + y (1) � u n +1 = . 3 This is an implicit Euler step (which is first order, but has nice stability properties) followed by a linear combination of some steps, which raises the order. The motivation for this approach is that we have legacy codes which have an implicit Euler scheme. Adding the ”filtering” step allows us to raise the order without touching the delicate workings of the implicit legacy code. Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 4 / 26

Generalization of Layton’s filtered Euler We can generalize Layton’s filtered Euler method by adding a ”pre-processing” as well as a ”post-processing” step: du n − 1 + (1 − d ) u n u n ˜ = u n + ∆ tF ( y (1) ) y (1) = ˜ 1 � − u n − 1 + 2(1 − d ) u n + 2 y (1) � u n +1 = . 3 − 2 d This extra freedom allows an entire family of second order methods that have implicit Euler as the driver. Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 5 / 26

Generalization of Layton’s filtered Euler With d = 0, this becomes Layton’s implicit Euler filter show above: u n + ∆ tF ( y (1) ) y (1) = y (1) − 1 � u n − 1 − 2 u n + y (1) � u n +1 = . 3 With d = 1, this becomes: u n − 1 + ∆ tF ( y (1) ) y (1) = 2 y (1) − u n − 1 u n +1 = Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 6 / 26

Numerical results To compare the effect of these filters, consider the problem dy dy = − 10 y 2 . IE + P 2 IE + P 2 ∆ t Implicit Euler Layton Filter LF Filter 5.0e-03 4.2642e-04 – 4.6646e-05 – 1.0503e-05 – 2.5e-03 2.1692e-04 0.97 1.1608e-05 2.00 2.7556e-06 1.93 1.25e-03 1.0944e-04 0.98 2.9001e-06 2.00 7.0661e-07 1.96 8.33e-04 7.3182e-05 0.99 1.2890e-06 1.99 3.1679e-07 1.97 6.25e-04 5.5000e-05 0.99 7.2557e-07 1.99 1.8138e-07 1.93 Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 7 / 26

Numerical Results Implicit Euler −3 Implicit Euler + Layton Filter Implicit Euler + LF Filter −3.5 −4 −4.5 −5 −5.5 −6 −6.5 −7 −3.4 −3.2 −3 −2.8 −2.6 −2.4 −2.2 Implicit Euler compared to Layton’s filter ( d = 0) and the leapfrog filter ( d = 1) Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 8 / 26

General Linear Methods Both the original method and the family of pre- and post-filtered methods can be written as a GLM. k k − 1 s d i ℓ u n − k + ℓ + ∆ t y n � � a i ℓ F ( u n − k + ℓ ) + ∆ t � a ij F ( y n i = ˆ j ) 1 ≤ i ≤ s ℓ =1 ℓ =1 j =1 k k − 1 s u n +1 = θ ℓ u n − k + ℓ + ∆ t � � b ℓ F ( u n − k + ℓ ) + ∆ t ˆ � b j F ( y n j ) , ℓ =1 ℓ =1 j =1 the u n − k + j are the previous steps, while the y n j are intermediate stages used to compute the next solution value u n +1 . Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 9 / 26

Layton’s Method as a General Linear Method Our first example is Layton’s filtered implicit Euler scheme : u n + ∆ tF ( y (1) ) y (1) = y (1) − 1 � u n − 1 − 2 u n + y (1) � u n +1 = . 3 We can write this as a one-stage two-step GLM: u n + ∆ tF ( y (1) ) y (1) = − 1 3 u n − 1 + 4 3 u n + 2 u n +1 3∆ tF ( y (1) ) , = which has the Butcher GLM tableau: � 2 − 1 4 � � � � � � � D = 0 1 , Θ = , A = 1 , b = 3 3 3 Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 10 / 26

Generalized Layton’s Method as a General Linear Method The family of second order methods with implicit Euler as the driver: du n − 1 + (1 − d ) u n u n ˜ = u n + ∆ tF ( y (1) ) y (1) = ˜ 1 � − u n − 1 + 2(1 − d ) u n + 2 y (1) � u n +1 = 3 − 2 d can be written as a GLM du n − 1 + (1 − d ) u n + ∆ tF ( y (1) ) y (1) = 2 d − 1 3 − 2 d u n − 1 + 4(1 − d ) 2 3 − 2 d u n + � y (1) � u n +1 = 3 − 2 d F , with tableau � � 4(1 − d ) 2 d − 1 2 � � � � � � D = d 1 − d , Θ = , A = 1 , b = 3 − 2 d 3 − 2 d 3 − 2 d Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 11 / 26

Generalized Layton’s Method as a General Linear Method The free parameter d gives us flexibility in the choice of filter. Two extreme cases to consider would be d = 0, which results in the Layton’s implicit Euler filter show above, and d = 1, which results in the scheme: u n − 1 + ∆ tF ( y (1) ) y (1) = 2 y (1) − u n − 1 u n +1 = which can be written in classic GLM form u n − 1 + ∆ tF ( y (1) ) y (1) = u n − 1 + 2∆ tF ( y (1) ) , u n +1 = with the Butcher tableau: � � � � � � � � D = 1 0 , Θ = 1 0 , A = 1 , b = 2 . Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 12 / 26

Producing Time Filters with a GLM approach The idea behind time filtering is that certain core methods are used, and are pre- and post-processed. Consider the core method with s − 1 stages and k steps: y n 1 = u n k − 1 k s d i ℓ u n − k + l + ∆ t � ˜ � � y n a i ℓ F ( u n − k + ℓ ) + ∆ t a ij F ( y n i = ˆ j ) , 2 ≤ i ≤ s ℓ =1 ℓ =1 j =1 k k − 1 s u n +1 = y n θ ℓ u n − k + ℓ + ∆ t ˜ � ˜ � ˆ � ˜ b ℓ F ( u n − k + ℓ ) + ∆ t b j F ( y n s = j ) j =1 ℓ =1 ℓ =1 Note that the final row coefficients are the same as the prior row coefficients d s ℓ , ˜ θ l = ˜ ˜ ˆ a s ℓ , ˜ b ℓ = ˆ b j = a sj . Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 13 / 26

Pre-filters using a GLM approach To pre-process the method, we modify u n . To do this, we change the initial stage so that u n is replaced by a linear combination of the old steps k 1 = u n − � d 1 ℓ u n − k + ℓ y n → y n 1 = ℓ =1 1 instead of u n in all the middle stages of the core method: and then use y n k − 1 k − 1 s d i ℓ u n − k + ℓ + ∆ t i = ˜ � ˜ � a i ℓ F ( u n − k + ℓ ) + ∆ t � y n d ik y n a ij F ( y n 1 + ˆ j ) ℓ =1 ℓ =1 j =1 k − 1 k s d i ℓ u n − k + l + ∆ t � � a il F ( u n − k + ℓ ) + ∆ t � a ij F ( y n = ˆ j ) , ℓ =1 ℓ =1 j =1 where d i ℓ = ˜ d i 1 d 1 ℓ + ˜ d i ℓ for ℓ = 1 , k − 1, and d ik = ˜ d i 1 d 1 k . Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 14 / 26

Post-filters using a GLM approach To post-process the method, we simply modify the final line to k k − 1 s u n +1 = θ ℓ u n − k + ℓ + ∆ t � � b ℓ F ( u n − k + ℓ ) + ∆ t ˆ � b j F ( y n j ) ℓ =1 ℓ =1 j =1 where the coefficients θ ℓ , ˆ b ℓ , b j can be chosen freely. 1 We can choose these coefficients so that the final stage can be written only as combinations of old steps u n − k + ℓ and stages y n j . The coefficients θ , ˆ b and b are referred to as the post-processing coefficients. 1 subject only to the order conditions Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 15 / 26

Time filtering using a GLM approach The overall method can now be written as k y n � d 1 ℓ u n − k + ℓ 1 = ℓ =1 k k − 1 s d i ℓ u n − k + ℓ + ∆ t � � � y n a il F ( u n − k + ℓ ) + ∆ t a ij F ( y n i = ˆ j ) 2 ≤ i ≤ s j =1 ℓ =1 ℓ =1 k k − 1 s u n +1 = θ ℓ u n − k + ℓ + ∆ t ˆ � � b ℓ F ( u n − k + ℓ ) + ∆ t � b j F ( y n j ) . ℓ =1 l =1 j =1 where d i ℓ = ˜ d i 1 d 1 ℓ + ˜ d i ℓ for ℓ = 1 , k − 1, and d ik = ˜ d i 1 d 1 k , and d 1 ℓ, θ ℓ , ˆ b ℓ , b j can be chosen for all ℓ and j . Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 16 / 26

What’s the point? Obviously, this is a fun game. But what does it give us? Approach the design of the pre- and post-processor as an optimization problem: Given the coefficients of the core method, we treat d 1 ℓ, θ ℓ , ˆ b ℓ , b j as free parameters The equality constraints are given by the order conditions we wish the new overall GLM to satisfy Additional requirements for low storage etc. can be imposed on the structure of the coefficients We define our objective function to be minimized (error constants) or maximized (linear stability region) We coded this optimization problem in MATLAB and quickly found methods that have desirable stability, accuracy, and efficiency properties. New GLMs based on trapezoid rule, BDF methods, fully implicit Runge–Kutta Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 17 / 26