We propose an algebraic basis for symmetric Strang splitting for first and second order accurate schemes for hyperbolic systems in N dimensions. Examples are given for two and three dimensions. Optimal stability is shown for symmetric systems. Lack of strong stability is shown for a non-symmetric example. Some numerical examples are presented for some Euler-like constant coefficient problems

Kucharik, Milan

Wendroff, Burton

UNT Digital Library

LA-UR- DJ-JJD/!J Approved for public release; distribution is unlimited. Tftle: Author(s): Intended for: MULTIDIMENSIONAL FIRST AND SECOND ORDER SYMMETRIC STRANG SPLITIING FOR HYPERBOLIC SYSTEMS MILAN KUCHARIK BURTON WENDROFF JOURNAL: APPLIED NUMERICAL MATHEMATICS A J LosAlamos NATIONAL LABORATORY --- EST.1943 --­Los Alamos National Laboratory, an affirmative action/equal opportunity employer, is operated by the Los Alamos National Security, LLC for the National Nuclear Security Administration of the U.S. Department of Energy under contract DE-AC52-06NA25396. By acceptance of this article, the publisher recognizes that the U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or to allow others to do so, for U.S. Government purposes. Los Alamos National Laboratory requests that the publisher identity this article as work performed under the auspices of the U.S. Department of Energy. Los Alamos National Laboratory strongly supports academic freedom and a researcher'S right to publish; as an institution, however, the Laboratory does not endorse the viewpoint of a publication or guarantee its technical correctness. Form 836 (7/06) Multidimensional First and Second Order Symmetric Strang Splitting for Hyperbolic Systems Milan Kucharika , Burton Wendroffb aTheoretical Division, T-5, Los Alamos National Laboratory MS-B284, Los Alamos, NM. 87545, USA, kucharik@lanl.gov, kucharik@karkulka.fjfi.cvut.cz b Theoretical Division, T-5, Los Alamos National Laboratory MS-B284, Los Alamos, NM. 87545, USA, bbw@lanl.gov Abstract We propose an algebraic basis for symmetric Strang splitting for first and second order accurate schemes for hyperbolic systems in N dimensions. Ex­amples are given for two and three dimensions. Optimal stability is shown for symmetric systems. Lack of strong stability is shown for a non-symmetric ex­ample. Some numerical examples are presented for some Euler-like constant coefficient problems. Key words: Multidimensional Hyperbolic Systems, Lax-Friedrichs, Lax-Wendroff, Symmetric Splitting, Optimal Stability 1. Introduction - One Dimension In his landmark paper [6] G. Strang introduced two splitting techniques designed to maintain second order accuracy for finite difference .methods for hyperbolic systems. The second of these has an advantage over the first in that it is symmetric with respect to interchange of dimension. In this note we shall develop symmetric splitting for both the first order Lax-Friedrichs type scheme, and for the second order Lax-Wendroff type, in N space dimensions, for hyperbolic conservation laws of the form N au _ L afj(u) at - j =1 ----a;;;-. Here, U = u(x,t), x = (Xl,··· ,XN), and fj(u) are vector valued func­tions. Preprint submitted to Applied Numerical Mathematics December 8, 2008 The difference schemes will be defined on uniform tensor product grids in terms of one-dimensional two point average and difference operators. In particular we define Qj == average over adjacent grid points in the j-th dimension and {3j == difference over adjacent grid points in the j-th dimension. The form of these operators in the more familiar grid index notation will be clear later. These operators all commute with each other. In these terms the familiar elementary one-dimensional finite difference operators are: Lax-Friedrichs (we suppress j here), and the two-step version of Lax-Wendroff, ~(h I + ApI 0 SF, where the A = Aj are the mesh ratios, The one dimensional time advance in the case of Lax-Wendroff is ",,,+1 _ ') Un u' - .. L , while in the case of Lax-Friedrichs it is SFUn n+1 S n+1 S S nU FU 2 = F 0 FU. 2 ,.. 2. 	 N-dimensional Lax-Friedrichs We offer two versions of N-dimensional Lax-Friedrichs type schemes. The first is ···0 (ap(N) + APt) (3P(NJiP(N)) ' (1) where P is the set of all permutations of {I ... N}. The second version is In either case, these operators map functions defined on the original grid to functions on the fully staggered grid at the half time step, or they map functions defined on the fully staggered grid to the original grid. So, as in one dimension, 2.1. 	 Linear Equivalence The two versions are linearly equivalent, that is, if the functions fj are replaced by multiplication by constant matrices then the two expressions are algebraically equivalent. To prove this, consider 1 A = N! 2:)ap(l) + Bp(l)) (ap(2) + Bp(2)) ... (ap(N) + Bp(N)), pEP where the B's are linear operators that commute with the a's. Then A has the form N N A = II a, + Ck 2:)L Bli ... B'k II al). (2) 1=1 1#1; 3 There are (~) ways to place the k factors containing k particular B's in the same order somewhere in an ordered product with the other factors, and there (N - k)! arrangements of the other factors each producing the same product of Q'S, so (Or, by symmetry, there must be N!jk! terms in the sum that are the product of k chosen B's in some order.) Next, collect all terms in (2) that have Bi as the leftmost factor, then all terms that have B2 as the leftmost, etc, to obtain which is the desired result. 2.2. Example In [4] the authors developed a 2 and 3 dimensional Lax-Friedrichs scheme, following the idea proposed in [2] of using iterated Riemann solvers. The second scheme above is the extension to ~ dimensions of what is found in [4]. Consider three dimensions. In order to keep the notation consistent with :4], let the system of differential equations be 8u 8g(u) 8h(u) at ax + 7iiJ + ---a;;-' vVith standard notation, first observe that the various products of Q'S are just the standard 2, 4 or 8 point averages. For example, 1 ( QiQ 3U ) = -(Ui 1" k + Ui+11"k + Ui 1" k+l + Ui+l 1" k+l)4 1,) , '.' '.j 4 Thus, 1 ( n8 	Ui,j,k + Ui+l,j,k + Ui,j,k+1 + Ui+1,j,k+l + Ui,j+1,k + Ui+l,j+l.k + Ui,j+l,k+l + Ui+1,j+l,k+1) + ~1 (f(U:+l,j+~,k+~) - f(U;,j+~'k+~)) + ~2 (g( V~+~,j+l,k+~) - g(v~+~,j,k+~) ) + ),3 (h(U/+1 '+1 k+J) h(W*+1 "+1 k)) .2 !2'}2" '2'}." For u* we have + ),2 (g(V**+l k+ 1 ) - g(v** k+ 1 ))4 2,)) 2 1,), 2 3+ ),4 (h(W**+l k+l) - h(W**+l k))~J 	 2' . IJ 2' and 1 ),32 (Hi,j,k+1 + Ui,j,k) + 6(h(Ui,),k+d - h(Ui,j,k)) 1 ),3 W" ** 1 k 2(Ui,j+l,k + Ui,j,k) + 6 (h(Ui,J+1,k) - h(Ui,j,k)).l'}+2' The functions v* and w* are similarly defined. 2.3. Stability If we consider now linear hyperbolic systems, that is N '" A. auL.. Jax j=l J Hyperbolicity means that every linear combination of the matrices Aj with real coefficients has only real eigenvalues. Strict hyperbolicity means that every linear combination of the matrices Aj with real coefficients has only real eigenvalues and is diagonalizable. \Ve are only able to establish a stability condition for symmetric hyperbolic systems 5 U sing the usual Fourier analysis, the symbol of the difference operator S is and ]V II,§II ~ I1llcos Oj1 + iAj sin OjAJ. j=l A sufficient stability condition then is that If the Aj are symmetric then Thus, for symmetric systems we have shown that the N-dimensional version of Lax-Friedrichs described above is stable if maxAj (spectral radius of AJ ~ 1. In other words, in this case the scheme is optimally stable since the CFL limit is 1. This is also true if each Aj is unitarily similar to a real diagonal matrix. 2.4. A counter-example Consider the linear system ap au ap-+­at ax ax au au at ax In this case the Lax-Friedrichs symbol is S cos 01 +iAA sin 0 6 ,.. where A (~ ~). Note that neither A nor!;' is diagonalizable, so that this is not a strictly hyperbolic system. Let us define z cosO + iAsinO, 0 < A S 1. Then, the n-th power of the Lax-Friedrichs symbol can be expressed as i n zn-l sin 0 ) !;'n ( o zn . The scheme is weakly stable, as guaranteed by the von Neumann condi­tion (see [5]), but not strongly stable. To see this, set and choose the sequence On such that . 20SIn n = 1-. n Then which is unbounded as n 00. 3. N-dimensional 2nd Order Scheme In order to develop an N-dimensional symmetric scheme of Lax-vVendroff type, as in [3], we need only replace the factors in (1) by 1+ A(3q where A q f(a + -(3f)2 7 or s ]]]. (4) 3.1. Accuracy S, as defined in (4) provides a second-order accurate difference scheme for the hyperbolic system. vVe give the basics of the proof. First, we can drop all terms after the factor i since they are O(ilt)a. Next, the central averaging operators can be replaced by the identity, and the central differences can be replaced by derivatives. since the errors made are of order (ilX)2ilt. Now from we have Let We need to show that for any smooth function u, But hfk (uk*+ ~t 8~k h(uk*)) = fk(U ) + ~t f~(u) 8~k h(u) + O(ilt)2, 8 and fk( llk*) fk( 1l) + ~t f~(1l) L j# from which (5) follows. 3.2. Example For example, in two dimensions the scheme is ll~J+l 1l7J· + Al (f(ll:+!) - f(ll*_!)) + A2 (g(V: ,+!) - g(v* ,_!)) ," I 2,J ' 2'] ',] 2 ',j 2 where 1 ( ** **) Al (f( ** )'2 lli+1,j + ll i ,j + 2 ll i+1,j and where 2 1 (lli,j+l + lli,j) + '21 A2 (g( lli,j+l g( lli,j)) . Similar expressions hold for v*. This agrees exactly with the 2D scheme in [3]. The 3D version obtained from (4) is not exactly the same as the one in [3]. 3.3. Stability Just as for Lax-Friedrichs, we can now establish the stability of the multi­dimensional scheme if the fj are symmetric matrices A j , for then for the symbol of (4) we have IISII IIIII + 2iAj sinOjAj(cosOj + iAsinOjAj)ll, where III + 2iA sinOA( cos 0 + iA sin OA)11 2 sup ((I-4A2A2sin20+4A4A4sin40+4A2A2 ocos2O)u, 1l) lIul!=l = 1 + (4A2 sin4 0) sup ((A2 A4 - A2)1l, 1l). (6) Ilull=l 9 So we have stability if (4) holds, that is, the scheme is optimally stable. REMARK: A simpler but non-optimally stable scheme is the following extension of the Richtmyer two-step scheme to N dimensions. It is shown in [7] that the CFL limit for this scheme is 1/VN. 4. Euler Equations We know that both 3D Lax-Friedrichs and the second order scheme are optimally stable for symmetric systems. Here we examine these features for the non-symmetric system of Euler equations. Our approach is based on numerical sampling of each scheme's amplification matrix, both as regards its eigenvalues and powers. We consider the 3D Euler equations in the form U t = f(U)x + g(U)y + h(U)z, where p pu pv pw E U= and pu pv pw pu2 + p puv puw f(U) = puv , g(U) = pv2 + p , h(U) = pvw puw pvw pw2 +p u (E + p) v (E + p) w (E + p) Here, p stands for fluid density, u, v, and ware components of velocity in each coordinate w = (u, v, w), E is a specific total energy defined by the 10 ideal gas equation of state + ~ P Ilw11 2 , p stands for pressure, and T is the gas adiabatic index. This system can be rewritten in a linearized form U t. = A U x + BUy + CUz , where A, B, and Care Jacobians of the flux functions, for example 8£(U)A 8U o o (1-1')'11 (1 1') W v o o U 1')/1,11 (l-1')uw and similarly for Band C. Note that the eigenvalues of this matrix are 1]" 1], + c, and 1], c, where c = JT p/p is the fluid sound speed. We have constructed the symbol (amplification matrix) of the 3D Lax­Friedrichs and the second order scheme. We only verify the von Neumann necessary stability condition, namely, that the maximum absolute value of the eigenvalues of the symbol is ::;: 1 for CFL::;: 1. If the maximal eigenvalue is greater than 1 the scheme is unstable. \Ve have chosen several values of interest for each quantity to be used in sampling. For each velocity component, we use all values from the following: 0, 1, ±2, ±3, and ±4, for sound speed we used values from: 0.1, 1, 2, 2.5, 5, and 10. It means that we sampled values for Mach number between 0 and 40 J3 ~ 69.3. As for CFL, we used values 1/10, 1/4, 1/2, 1-1/10, 1 1/100, 1 -10-5 , 1, 1+ 10-5, 1 + 1/100, 1 + 1/10, and 2. Fourier angles were chosen to cover all quarters and thirds of 7r in the whole 27r range, and monatomic gas constant ~I = 5/3 was used. The maximal eigenvalue for each scheme as a function of CFL is presented in Table 1. For computation of eigenvalues the LAPACK [1] package was used. This verifies, at least for a supposedly typical set of arguments of the sym­bol, that Lax-Friedrichs and our second order scheme are optimally (weakly) stable even for non-symmetric systems. 4.1. Zem Sound Speed Included in the results of Table 1 the smallest value of the sound speed is c .1. We consider here what happens if we include c O. 11 CFL Lax-Friedrichs 2nd order scheme 1 1/10 1.0000000000 1.0000000000 1 1/100 1.0000000000 1.0000000000 1 10-5 1.0000000000 1.0000000000 1 1.0000000000 1.0000000000 1 + 10-5 1.0000100000 1.0001200054 1 + 1/100 1.0100000000 1.1255130848 1 + 1 1.2898723387 2.8632880000 Table 1: Maximal of amplification matrix for 3D Lax-Friedrichs and the second order scheme for several CFL numbers. First, for the same sampling set, the schemes are optimally weakly stable. Second, for CFL=l, the maximum component of the n-th power of the symbol samples increases linearly with n. This was tested for powers n (2,4,· .. ,215 ). The maximum is taken over the hydrodynamic data samples and Fourier angles. Third, following section 2.4, for CFL< 1 and sin2 OJ = l/n, we again observe that the powers of the symbol are not bounded, Fig. 1. Acknowledgments This work was performed under the auspices of the :.Jational :.Juclear Security Administration of the US Department of Energy at Los Alamos l\ational Labora­tory, under Contract W-7405-ENG-36 and Contract DE-AC52-06NA25396. The authors acknowledge the partial support of the DOE Advance Simulation and Computing (ASC) Program and the DOE Office of Science ASCR Program, and the Laboratory Directed Research and Development program (LDRD) at the Los Alamos National Laboratory. References [1] Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users' Guide, 3rd ed., Society for Industrial and Applied Math­ematics, Philadelphia, PA, 1999. [2] 	 T. Boukadida, A.-y' LeRoux, A new version of the two-dimensional Lax­Friedrichs scheme, Math. Compo 35 (1994) 541-553. 12 3D LF scheme 7 6 5 ~:: 4 x III .s 03 rj 7 6 -5 c::: sa. 4 ~ .s ~ 3 '" g 2 6 8 10 12 14 I09 (n) 23D 2nd order scheme - 0.1 - 0.25 0.45 - 0.5 - 0.55 - 0.9 - 0.99 - 0.99999 - 1 - 1.00001 -1.01 - 1.1 - 2 16 18 - 0.1 - 0.25 0.45 - 0.5 - 0.55 - 0.9 - 0.99 - 0.99999 - 1 - 1.00001 - 1.01 - 1.1 - 22 I 12 14 16 182 4 6 Figure 1: Growth rate (maximal component of symbol sn) for 3D first order Lax-Friedrichs and second order Lax-Wendroff schemes for different CFL numbers and powers n . 13 [3] 	 YI. Kucharik, R. Liska, S. Steinberg, B. Wendroff, Optimally-stable second~order accurate difference schemes for non-linear conservation laws in 3D, Applied Numerical Mathematics 56 (2006) 589 - 607. [4] 	 R. Liska, Wendroff, Composite sehemes for eonservation laws, SIAM J. Numer. Anal. 35 (1998) 2250 2271. [5] 	 R. D. Richtmyer, K. \V..~vlorton, Difference Methods for Initial Value Problems, John Wiley, 1967, 2nd Ed. [6: 	 G. Strang, On the construction and comparison of difference schemes, SIAM J. Numer. Anal. 5 (1968) 506 517. [7] 	 B. Wendroff, Well-posed problems and stable difference operators, SIAM J. Numer. Anal. 5 (1) (1968) 71-82. 14 

Multidimensional first and second order symmetric strang splitting for hyperbolic systems

https://digital.library.unt.edu/ark:/67531/metadc935337/m2/1/high_res_d/956668.pdf

Multidimensional first and second order symmetric strang splitting for hyperbolic systems

Abstract

Similar works

Full text

Available Versions

UNT Digital Library