Positivity-Preserving High-Order Discontinuous Galerkin Methods: Implicit Time Stepping and Applications to Relativistic Hydrodynamics by Tong Qin B.Sc., University of Science and Technology of China, Hefei, P. R. China, 2011 M.Sc., University of British Columbia, Vancouver, BC, 2013 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2017 c Copyright 2017 by Tong Qin This dissertation by Tong Qin is accepted in its present form by the Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Chi-Wang Shu, Ph.D., Advisor Recommended to the Graduate Council Date Johnny Guzm´an, Ph.D., Reader Date Jerome Darbon, Ph.D., Reader Approved by the Graduate Council Date Andrew G. Campbell, Dean of the Graduate School iii Vita Education Brown University, Providence, RI, U.S.A., Ph. D. in Applied Mathematics, May, 2017 (expected), Advisor: Prof. Chi-Wang Shu. University of British Columbia (UBC), Vancouver, Canada, M. Sc. in Applied Mathematics, Nov., 2013, Advisor: Prof. Dominik Sch¨otzau. University of Science and Technology of China (USTC), Hefei, China, B. Sc. in Mathematics and Applied Mathematics, Jul. 2011, Advisor: Prof. Yan Xu. Publications and Preprints 1. T. Qin and C.-W. Shu, Implicit positivity preserving high order discontinuous Galerkin methods for conservation laws, submitted. 2. T. Qin, C.-W. Shu and Y. Yang, Bound-preserving discontinuous Galerkin methods for relativistic hydrodynamics, Journal of Computational Physics, vol. 315, pp. 323-347, 2016. 3. R. Oyarz´ua, T. Qin and D. Sch¨otzau, An exactly divergence-free finite ele- ment method for a generalized Boussinesq problem, IMA Journal of Numerical Analysis, vol. 34 (3), pp. 1104-1135, 2014. iv Honors and Awards • 2010, Outstanding Undergraduate Student Research Project, USTC. • 2007-2011, Outstanding Student Scholarship, USTC. • 2010, Second prize in Statistical Modelling Competition, USTC. • 2007, Second prize in National Olympia Physics Competition. Teaching Experience At Brown • Spring 2015, Teaching Assistant, APMA 1200 Operational Analysis: Probabilistic Models. • Fall 2014, Teaching Assistant, APMA 0360 Methods of Applied Mathematics. At UBC • Spring 2013/2012, Teaching Assistant, MATH 210 Introduction to Mathematical Computing. • Fall 2012/2011, Teaching Assistant, MATH 255 Ordinary Differential Equations. • 2011-2013, Graduate Tutor in the Math Learning Center. v Acknowledgments First of all, I would like to express my deepest gratitude to my supervisor Prof. Chi-Wang Shu for his advice and encouragement when I got stuck and depressed, for his bluntly pointing out of my shortcomings, for his generous sharing his life stories and experience with me and for his thoughtful advice for both of my career and my family life. He not only showed me how to do high-standard and useful research but also taught me how to work efficiently and how to be a self-disciplined, responsible and integral person. The most valuable thing I learned from him is to be generous to others while being strict with ourselves. All in all, he has set me a life-long role model. Secondly, I would like to thank Prof. Johnny Guzm´an and Prof. Jerome Dar- bon who have graciously agreed to spend their time in reading this dissertation and providing invaluable feedbacks. I also owe my thanks to all the other professors who have once shared their wisdoms with me, either in class or in research projects, in- cluding Prof. Mark Ainsworth, Prof. Yan Guo, Prof. Walter Strauss, Prof. Hongjie Dong, Prof. Paul Dupuis at Brown as well as Prof. Dominik Sch¨otzau, Prof. Michael Ward and Prof. Uri Ascher at UBC, Prof. Zhiming Chen at CAS and Prof. Yan Xu at USTC. Thirdly, my thanks also go to Dr. Xiangxiong Zhang and Dr. Yang Yang. Dr. Xiangxiong Zhang’s intellectual work in his thesis laid a foundation for this vi dissertation. The second part of the dissertation is based on the work [55], in which Dr. Yang Yang has made a great contribution. Moreover, their stories at Brown constantly remind me how to be a hard-working and resilient student. Next, all my friends at Brown have helped me a lot during my four years’ grad- uate study. Among them I would thank Chengke Shi, for numerous enthusiastic discussions on math, life and many other topics. Thanks to Zheng Sun, for being a mirror reflecting many of my shortcomings. Thanks to Tianheng Chen, Guosheng Fu and Juntao Huang, for sharing their research with me and for their help with my life. Thanks to Weifeng Sun for being a considerate roommate and a reliable friend. Thanks to Bingyu Zhao and Kunrui Wang, for organizing our weekly academic fam- ily gathering. Thanks to my academic buddy Christian Glusa, for his advice on my preliminary exam and for many conversations on each other’s research. Thanks to Dr. Jeffery Miller for kindly sharing this thesis template around. At last, thanks to Prof. Henry Majewski for hosting me and providing me with help during my very first year in Providence. All these could not happen without the generous support from the Brown Univer- sity and the Graduate School at Brown as well as the service of all the staff members of DAM at Brown. Last but not the least, I will keep my special gratitude to my family. Thanks to my parents for their unconditional support and sacrifice. Thanks to my wife, my beautiful sunshine, Mrs. Yang Chen, for her constant support and tolerance. vii Abstract of “Positivity-Preserving High-Order Discontinuous Galerkin Methods: Im- plicit Time Stepping and Applications to Relativistic Hydrodynamics”, by Tong Qin, Ph.D., Brown University, May 2017 The positivity-preserving property is a highly desirable property when designing high order numerical methods for hyperbolic conservation laws, since negative val- ues sometimes cause ill-posedness of the problem and blow-ups of the algorithms. The general framework for constructing positivity-preserving schemes for solving hy- perbolic conservation laws have been proposed in (X. Zhang and C.-W. Shu, Journal of Computational Physics, 229 (2010), pp. 3091–3120) and (X. Zhang and C.-W. Shu, Journal of Computational Physics, 229 (2010), pp. 8918–8934). In this disser- tation, we extend this framework to DG methods with implicit discretizations and to DG methods for solving the relativistic hydrodynamics (RHD). Due to the the Courant-Friedrichs-Levis (CFL) number constraint, explicit DG methods are impractical for problems involving unstructured and extremely varying meshes or long-time simulations. Instead, implicit DG schemes are often popular in practice, especially in the computational fluid dynamics (CFD) community. In the first part of this dissertation, we develop a high-order positivity-preserving DG method with the backward Euler time discretization for solving conservation laws, basing on a generalization of the Zhang-Shu positivity-preserving limiter. Both the analysis and numerical experiments indicate that a lower bound for the CFL number is required to obtain the positivity-preserving property for the numerical schemes. The proposed method not only preserves the positivity of the numerical approxima- tion without compromising the designed high-order accuracy, but also helps accel- erate the convergence towards the steady-state solution and add robustness to the nonlinear solver. For RHD systems, the density and the pressure are positive physical quantities and the velocity is bounded by the speed of light. The violation of these bounds will result in ill-posedness of the problem and blow-up of the code, especially in extreme relativistic cases. It is usually hard to maintain these physical bounds without sac- rificing the numerical accuracy. In the second part, we develop a bound-preserving DG method to solve RHD systems by extending the bound-preserving limiter for the non-relativistic hydrodynamics. The proposed method has the following features. It can theoretically guarantee to preserve the physical bounds for the numerical ap- proximation and maintain its designed high order accuracy. Moreover, it renders L1 -stability to the numerical scheme. The robustness of the scheme is tested on various extreme relativistic cases, including relativistic jets. Contents Vita iv Acknowledgments vi 1 Introduction 1 2 DG Methods and Framework for Bound-Preserving Techniques 5 2.1 DG methods for conservation laws . . . . . . . . . . . . . . . . . . . . 6 2.1.1 Spacial discretizations . . . . . . . . . . . . . . . . . . . . . . 7 2.1.2 Time discretizations . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 General framework for positivity and bound preserving techniques . . 9 3 Positivity-Preserving DG Methods with Implicit Time Stepping for Conservation Laws 12 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 Positivity-preserving backward Euler DG scheme for scalar equations 16 3.2.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.2 CFL condition for P 1 -DG methods on uniform meshes . . . . 20 3.2.3 CFL condition for arbitrary polynomial degree on general meshes 26 3.2.4 Scaling limiter . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.5 Algorithm for scalar equations . . . . . . . . . . . . . . . . . . 42 ix 3.3 Positivity-preserving backward Euler DG scheme for compressible Eu- ler systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.4 Positivity-preserving Crank-Nicolson DG scheme . . . . . . . . . . . . 46 3.4.1 CFL condition for P 1 -DG on uniform meshes . . . . . . . . . 47 3.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5.1 Accuracy tests . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.5.2 Moving shocks . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.5.3 Compressible Euler system . . . . . . . . . . . . . . . . . . . . 66 3.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4 Positivity and Bound-Preserving DG Methods for Relativistic Hydrodynamics 71 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.2 Numerical algorithm in one space dimension . . . . . . . . . . . . . . 78 4.2.1 The DG scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 78 4.2.2 Bound-preserving technique . . . . . . . . . . . . . . . . . . . 79 4.2.3 High order time discretizations . . . . . . . . . . . . . . . . . 90 4.3 Numerical algorithm in two space dimensions . . . . . . . . . . . . . 91 4.4 Application to relativistic jets . . . . . . . . . . . . . . . . . . . . . . 98 4.5 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 102 4.5.1 One-dimensional experiments . . . . . . . . . . . . . . . . . . 103 4.5.2 Two-dimensional experiments . . . . . . . . . . . . . . . . . . 111 4.6 Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5 Conclusion 120 x List of Tables 3.1 Values of rk for k = 1, · · · , 5 and for Legendre-Gauss-Lobatto (LGL) and Legendre-Gauss quadrature rules respectively. . . . . . . . . . . . 40 3.2 Minimum cell average after one time step for odd k. . . . . . . . . . . 40 3.3 Minimum cell average after one time step for even k. . . . . . . . . . 41 3.4 Values of minj u¯1j with initial condition (3.45) for N = 4, 8, 16, 32 and different λ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.5 Values of minj u¯1j with initial condition (3.46) for N = 4, 8, 16, 32 and different λ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.6 Error table for Example 3.5.1.1 at T = 1 with λ = 0.81, uniform mesh. 58 3.7 Error table for Example 3.5.1.2, approximation of the steady state solution to the linear problem (3.48), without the positivity preserving limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.8 Error table for Example 3.5.1.2, approximation of the steady state solution to the linear problem (3.48) with the positivity preserving limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.9 Error table for Example 3.5.1.3, approximation of the steady state solution to the Burgers equation (3.49), without the positivity pre- serving limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 3.10 Error table for Example 3.5.1.3, approximation of the steady state solution to the Burgers equation (3.49), with the positivity preserving limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 xi 3.11 Error table for Example 3.5.1.3, approximation of the steady-state so- lution to the Burgers equation (3.50), with and without the positivity preserving limiter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.1 Example 4.5.1.1: One-dimensional accuracy test at T = 0.4 for the second-, third- and fourth-order DG methods with and without the BP limiters. The CFL number for the multistep method is one-third of that for the RK method. k is the degree of polynomial and h is the meshsize. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 4.2 Example 4.5.1.1: Two-dimensional accuracy test at T = 0.2 for the third-order RKDG method with the BP limiter. k is the degree of polynomial and h is the mesh size. . . . . . . . . . . . . . . . . . . . 114 xii List of Figures 3.1 Example 3.5.2.1: at T = 0.2, with ∆t = 0.266 maxj hj and N = 120. Solid line is the exact solution. Blue squares are point values on the Legendre-Gauss-Lobatto points. . . . . . . . . . . . . . . . . . . . . . 64 3.2 Example 3.5.2.2 at T = 1.5. Solid line is the exact solution. Blue squares are point values of the numerical approximation on the Legendre-Gauss-Lobatto points. Mesh with N = 120 cells and CFL = 2. 65 3.3 Example 3.5.2.3 at T = 0.6. Solid line is the exact solution. Blue squares are point values of the numerical approximation on the Legendre-Gauss-Lobatto points. Mesh with N = 120 cells and CFL = 3.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.4 Example 3.5.3.1 at T = 0.01 with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are numerical approximations. . . . . . . . . . . . . . . . 67 3.5 Example 3.5.3.2 at T = 0.1, with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are the numerical approximations. . . . . . . . . . . . . 68 3.6 Example 3.5.3.3 at T = 0.003, with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are the numerical approximations. . . . . . . . . . . . . 69 xiii 4.1 Plots for the lower bound of α in (4.19) for the bound-preserving requirement (denoted by αbp ) and the spectral radius (4.28) of the Jacobian matrix (denoted by αstandard ). The left panel is to keep γ = 20 constant and the right one is to keep h = 20 constant. . . . . . 85 4.2 Example 4.5.1.1: Log-log plot of the CPU time against the L2 error for polynomial degrees k = 1, 2, 3 with and without the BP limiter. . 105 4.3 Example 4.5.1.2: Blast wave with initial condition (4.51) at T = 0.5 on the mesh of 200 cells, approximated by the third-order RKDG method with the BP limiter. The squares represent the approximate cell averages and the solid line is the exact solution. . . . . . . . . . . 106 4.4 Example 4.5.1.3: Strong blast wave with initial condition (4.52) at T = 0.4 approximated by the third-order RKDG method with the BP limiter on the mesh with 200 cells. The square represents the approximate cell averages and the solid line is the exact solution. . . . 107 4.5 Example 4.5.1.3: Strong blast wave with initial condition (4.53) at T = 0.3 approximated by the third-order RKDG method with the BP limiter. The mesh is decomposed into 400 cells (left column) and 800 cells (right column). The square represents the approximate cell averages and the solid line is the exact solution. . . . . . . . . . . . . 108 4.6 Example 4.5.1.4: One-dimensional shock-heating problem on the mesh of 200 cells, at T = 1.5 with vin = 0.999999, approximated by the third-order RKDG method with the BP limiter. The squares represent the numerical results and the solid line is the exact solution. . . . . . 109 4.7 Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity at T = 0.6, approximated by the third-order RKDG method with the BP limiter. Approximations of the den- sity and the transverse velocity on meshes of 400 and 6400 cells. The squares represent the approximate cell averages and the solid lines are the exact solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 xiv 4.8 Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity at T = 0.6, approximated by the third-order RKDG method with the BP limiter. The approximations of w and p on meshes of 400 and 6400 cells are presented respectively. The squares represent the approximate cell averages and the solid lines are the exact solutions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.9 Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity v = 0.999 at T = 0.6, approximated by the third- order RKDG method with the BP limiter. The approximations of the w, v, p and ρ on 3200 cells are presented. The blue solid lines are the numerical results and the black ones are the exact solutions. . . . . . 113 4.10 Example 4.5.2.2: Blast wave with initial condition (4.56) propagating √ 2 2 in [0, 2 ]. At time T = 0.3 on a 120 × 120 mesh. In the right column, solid lines are exact solution and the blue squares are the approximate ones. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 4.11 Example 4.5.2.3: 400 × 400 cells at T = 0.7 approximated by the third-order RKDG method with the BP limiter. Thirty equally spaced contours of the logarithm of the proper density are plotted. . . . . . 116 4.12 Example 4.57: Relativistic jets of model C2 with initial condition (4.57). The left panel is at T = 20 and the right is at T = 35. The resolution is 10 points per jet radius. . . . . . . . . . . . . . . . . . . 117 4.13 Example 4.5.2.4: Relativistic jet of the model C3 with initial condition (4.58) at T = 30, approximated by RKDG with the BP limiter. The resolution is 10 points per jet radius. . . . . . . . . . . . . . . . . . . 118 xv Chapter One Introduction 2 The discontinuous Galerkin (DG) methods were first introduced by Reed and Hill [59] for solving neutron transport equations on triangular meshes. Afterwards, Cockburn et al. developed the Runge-Kutta DG (RKDG) methods in a sequence of papers [12, 11, 10, 13] for solving the hyperbolic conservation laws. After that, DG methods have seen great success in many areas of applications, including the computational fuild dynamics (CFD), semi-conductor device design, computational comsmology and so on. The DG method features mathematically provable high-order accuracy and stability, extremely local data structure, great parallel efficiency and flexibility for dealing with unstructured meshes and h-p adaptivity. One fundamental challenge in desgining high-order DG methods for hyperblic conservation laws is to maintain certain stability without compromising the designed high-order accuracy. This is usually done by limiting the numerical approximation or by modifying the schemes to let the numerical solution respect certain properties of the physical solution, such as the total variation diminishing or bounded (TVD/TVB) properties [12], entropy inequalities [8], maximum principle and so on. One of the recent great achievements in overcoming this difficulty is the genuinely high-order maximum-principle satisfying and positivity preserving DG methods developed by Zhang and Shu in [86] and [87]. In this dissertation, we generalize this work in the following two directions. 1. Generalize the positivity-preserving limiter for the explicit methods to the DG methods with backward Euler time discretizations [54]. 2. Generalize the bound-preserving limiter for non-relativistic hydoynmaics to the relativistic one [55]. The first topic is motivated by the practical application of DG methods in the computational fluid dynmaics (CFD) community. When solving the aerospace prob- lems, unstructured and extremely varying meshes are always adopted to resolve 3 complex geometries or boundary layers and long-time simulations are always needed for the steady-state simulation. The CFL constraint for explicit time discretizations make explicit methods, such as RKDG methods, impractical in these real applictions. To circumvent this difficulty, implicit methods, which have much larger stability re- gions, are often popular in practice. In particular, for DG methods, in [29] Jiang and Shu have shown by the cell entropy inequality that a class of implicit time dicretiza- tions, including the backward Euler ant the Crank-Nicolson methods, are uncondi- tionally L2 stable. However, only a few works can be found in the literature dealing with the positivity-preserving property for high-order DG methods coupled with im- plicit time discretizations. In the first part of this dissertation, we introduce how to construct positivity-preserving DG methods for conservation laws with arbitrary order of spacial accuracy. Since our focus is on the positivity-preserving property of the spacial discretization, the time is discretized with the simple backward Euler methods. The Crank-Nicolson time discretization is also briefly considered for linear element on uniform meshes. High-order implicit time discretization will be explored in future work. We derive a CFL condition for the scalar linear equation under which the cell average is promised to be positive in each time step. Then the positivity of the whole polynomial, at least values at quadrature points, can be obtained by the standard scaling limiter [33, 86]. The second topic is to generlize the bound-preserving limiter for the compress- ible Euler system with and without source terms in [87, 89] to the RHD system in both of the Cartesian and cylindrical coordinates. One of the challenges in solving relativistic hydrodynamics is how to preserve the physical bounds of the solution. Compared with the non-relativistic one, for the relativistic hydrodynamics, the prim- itive variables and the conservative variables are related in a nonlinear way due to the Lorentz transformation. In order to recover the positive pressure from the conserva- 4 tive variables in each time step, a nonlinear equation needs to be solved. However, the existence of the positive root of this nonlinear equation requires the numerical approximation for the conservative variables to fall into certain admissible set. We develop a high-order bound-preserving DG method to solve this robutness problem in both one and two dimensions with both Cartisian and cylindrical coordinates respec- tively. The method is tested on extreme relativistic cases including the relativistic jets. This dissertation consists of the following chapters. In Chapter 2, we will do a brief review of the DG methods and the general idea of the positivity and bound- preserving techniques. Then Chapter 3 includes the extension of the positivity- preserving limiter to the implicit DG methods. The extension to the RHD system is introduced in Chapter 4. Concluding remarks will be given in Chapter 5. Chapter Two DG Methods and Framework for Bound-Preserving Techniques 6 In this chapter, we use the following one-dimensional conservation law system ut + f (u)x = 0, (x, t) ∈ Ω × [0, ∞), (2.1) u(x, 0) = u0 (x), x ∈ Ω, with appropriate boundary conditions to illustrate the formulation of the DG method and the general framework for the positivity and bound-preserving technique. Here Ω = [0, 2π], u : [0, 2π] × R+ → Rd is a vector of state variables and the function f : Rd → Rd is the flux function. When d = 1, we use the non-bold face letters u and f to denote the scalar variable and the flux function. The well posedness of the problem (2.1) requires that for each fixed u, the Jaco- ∂f bian matrix ∂u has real eigenvalues {ζj (u)}dj=1 , and has a complete set of d linearly independent eigenvectors. 2.1 DG methods for conservation laws To define the DG scheme for (2.1), let us first fix some notations. We decompose the domain Ω = [0, 2π] into N subintervals, Ij = [xj− 1 , xj+ 1 ], for j = 1, 2, · · · , N . The 2 2 size of each subinterval is denoted by hj . Define Iˆ = [−1, 1] to be the reference cell and define Tj (x) = 2(x − xj )/hj to be the affine mapping between the intervals Ij ˆ where xj = (x 1 + x 1 )/2 is the midpoint of Ij . Moreover, let (·, ·)j denote and I, j− j+2 2 the usual L2 inner product on Ij , (·, ·)Iˆ the one on Iˆ and let (·, ·)Ω denote the inner product on Ω. Then we define the approximation space Vh = {v ∈ L2 (Ω) : v|Ij ∈ Pk (Ij ), ∀j = 1, · · · N } 7 where Pk (Ij ) denotes the polynomial space on Ij with degree up to k and its vector version Vh = [Vh ]d . 2.1.1 Spacial discretizations For fixed t, the spacial discretization with the DG method for (2.1) is defined by seeking the approximation uh (t) ∈ Vh , such that in each subinterval Ij , we have d (uh (t), v)j − (f (uh (t)), vx )j + ˆfj+ 1 (uh (t))v(x− j+ ˆ 1 ) − fj− 1 (uh (t))v(x + j− 12 ) = 0, (2.2) dt 2 2 2 ∀v ∈ Vh , where ˆfj+ 1 (u) = ˆf (u(x− j+ 1 ), u(x+ j+ 1 )) and ˆf (·, ·) is an exact or approximate 2 2 2 Riemann solver [69] for system and a monotone flux [32] for scalar equation. There are many choices for ˆf (see [69, 32]), for example, the global Lax-Friedrichs flux in the following form ˆf (a, b) = 1 [f (a) + f (b) − α(b − a)], (2.3) 2 where α = maxi kζi (u(·, t))k∞ . 2.1.2 Time discretizations The semi-discrete scheme (2.2) can be written in the following compact form, d (uh (t), v)j = Lj (uh (t), v), ∀v ∈ Vh , ∀j = 1, . . . , N (2.4) dt with the shorthand notation h i Lj (uh (t), v) = (f (uh (t)), vx )j − ˆfj+ 1 (uh (t))v(x− j+ 1 ) − ˆ fj− 1 (u h (t))v(x + j− 1 ) . (2.5) 2 2 2 2 8 For the ODE system (2.4), one can either use explicit methods or implicit ones to further discretize it. For the explicit methods, a common choice is the class of strong stability pre- serving (SSP) methods [64, 65, 23]. Among them the most popular ones are the following third-order SSP Runge-Kutta method [65], (1) (uh , v)j = (unh , v)j + ∆tLj (unh , v)j , ∀j, ∀v ∈ Vh , 3 1 (1) 1 (1) (u(2) , v)j = (un , v)j + (uh , v)j + ∆tLj (uh , v), ∀j, ∀v ∈ Vh , (2.6) 4 4 4 1 n 2 (2) 2 (2) (un+1 h , v)j = (uh , v)j + (uh , v)j + ∆tLj (uh , v), ∀j, ∀v ∈ Vh , 3 3 3 and the following third order SSP multistep method [64], 16 (un+1 h , v)j = [(un , v)j + 3∆tLj (unh , v)] 27 h  (2.7) 11 n−3 12 n−3 + (uh , v)j + ∆tLj (uh , v) , ∀j, ∀v ∈ Vh . 27 11 where unh ∈ Vh denotes the numerical approximation at time level n. For the implicit discretizations, many methods have been considered in practice, including the backward difference formula (BDF) methods [27], the implicit modified extended BDF (MEBDF) method [46], linearly implicit Rosenbrock-type RK method [2] and the fully implicit RK methods [28, 49]. In this dissertation, we considered the simplest one, i.e., the backward Euler method (un+1 n+1 n h , v)j − ∆tLj (uh , v) = (uh , v)j (2.8) 9 for all v ∈ Vh and j = 1, · · · , N . And the Crank-Nicolson method ∆t ∆t (un+1 h , v)j − Lj (un+1 n h , v) = (uh , v)j + Lj (unh , v) (2.9) 2 2 for all v ∈ Vh and j = 1, · · · , N . In the following, we use unj to denote unh |Ij , (unj+ 1 )± to denote unj (x± j+ 1 ) and use 2 2 ¯ nj u to denote the cell average of unj in the interval Ij . 2.2 General framework for positivity and bound preserving techniques The physical admissible solution to (2.1) usually falls into some admissible set G. In most applications, the set G can be verified to be a convex set. For example, for the scalar equation, it is well known that the entropy solution satisfies the following maximum principle min u0 (x) ≤ u(x, t) ≤ max u0 (x), ∀t ≥ 0. (2.10) x∈Ω x∈Ω In this case the set G is [m, M ] with m = minx∈Ω u0 (x) and M = maxx∈Ω u0 (x). Our goal is to design a high-order DG scheme such that unh (x) ∈ G, ∀x ∈ Ω =⇒ un+1 h (x) ∈ G, ∀x ∈ Ω (2.11) A general framework for constructing such DG schemes has been proposed in [86, 87], which can be summarized in two steps as below. 10 Step 1 First, in each cell Ij , given the approximation at time level n is admissible, i.e., unj (x) ∈ G for any x ∈ Ij , find a sufficient condition such that we have part of the approximation, i.e., the cell average at the next time level is in G, ¯ n+1 u j ∈ G. Step 2 Next, in each cell Ij , we limit the whole approximation un+1 j (x) into the set G by invoking the following scaling limiter [33, 86, 87], ˜ n+1 u j = θj [un+1 j ¯ n+1 −u j ¯ n+1 ]+u j (2.12) ˜ n+1 where θj ∈ (0, 1) is solved such that u j ∈ G. The main difficulty lies in the first step. The sufficient conditions for DG methods with explicit SSP time discretizations, have been derived in [86, 87, 89] for general scalar conservation laws and the compressible Euler system with and without source terms respectively. In Chapter 4, we will extend this procedure to the RHD system and design a high-order bound preserving explicit DG methods. For the implicit time discretization, which is popular in practical application, the positivity-preserving DG methods been less explored. In Chapter 3, we mainly investigate the simplest case, i.e., the backward Euler time discretization for the one dimensional conservation laws. And we also briefly consider the Crank-Nicolson time discretization for linear DG methods on uniform meshes. To achieve the goal in Step 1, a set of sufficient conditions are derived for the scalar linear equations. These conditions are further verified by numerical examples for nonlinear scalar equations. The applicability of the proposed method is also tested for the compressible Euler systems. As for the second step, the main issue is whether this modification will destroy 11 the original high-order accuracy of the DG method or not. For scalar equations, the complete result for arbitrary order of accuracy is first given in [86] and the proof is further refined in [85]. The result is included in the following lemma. Lemma 2.2.1. For scalar equations where G = [m, M ], the modified polynomial u˜n+1 j (x) as defined in (2.12) preserves the original order of accuracy in the following sense, un+1 |˜ j (x) − un+1 j (x)| ≤ Ck max |un+1 j (x) − u(x)| x∈Ij where u is the smooth solution to (2.1) with d = 1. Basically, what this lemma says is that for the scalar equation the error we commit in the limiting procedure is bounded by the error of the original approximation up to a constant depending on k. For the system, even though there is no rigorous result as Lemma 2.2.1, both truncation error analysis [88] and numerical results show that the scaling limiter (2.12) preserves the high-order accuracy as well. Chapter Three Positivity-Preserving DG Methods with Implicit Time Stepping for Conservation Laws 13 In this chapter, we consider the following conservation law ut + f (u)x = 0, (x, t) ∈ [0, 2π] × [0, +∞) (3.1) u(x, 0) = u0 (x), x ∈ [0, 2π] and its system version with appropriate boundary conditions. We focus on this one-dimensional case, even though the result can be easily generalized to multidi- mensional tensor product meshes and polynomial spaces. 3.1 Introduction For scalar conservation laws, it is well known that the entropy solution satisfies the following maximum principle min u0 (x) ≤ u(x, t) ≤ max u0 (x), ∀t ≥ 0. x∈[0,2π] x∈[0,2π] In particular, if the initial condition is positive, then the entropy solution should satisfy the following positivity-preserving property u0 (x) ≥ 0 =⇒ u(x, t) ≥ 0, ∀t ≥ 0 (3.2) For systems, even though the entropy solution does not satisfy the maximum princi- ple in general, the physically relevant solution, for example the density and pressure in the compressible Euler system, is always positive. In this chapter, the words “posi- tive” and “positivity” are used actually to mean “nonnegative” and “nonnegativity”. We shall use “strictly positive” to mean the usual “positive”. When designing numerical methods, we would like our numerical approximations 14 to respect this positivity-preserving property (3.2), not only because it makes the nu- merical approximation physically meaningful, but also it makes the numerical scheme more robust, since negative values sometimes cause ill-posedness of the problem and blow-ups of the numerical algorithm [24]. In recent years, the positivity-preserving DG schemes have been actively designed and applied for solving hyperbolic conserva- tion laws [86, 87, 78, 79, 55, 9, 85]. All these methods are coupled with explicit tem- poral discretizations, such as strong stability preserving (SSP) Runge-Kutta (RK) methods [65, 23] and multi-step methods [64]. Explicit temporal discretizations enjoy many advantages, for example, the easiness in handling the nonlinear terms and boundary conditions, high-order accuracy with SSP properties [23], low stor- age requirement and so on. However, they suffer from the CFL constraint. For DG methods, to obtain the linear stability [1] or the maximum-principle stability [86], the CFL constraint becomes more and more severe as we increase the polynomial degree in the approximation space. Such stringent time stepping restriction makes explicit methods impractical in computations involving unstructured and extremely varying meshes, viscous effect [51] or long-time simulations for steady-state calculation [27]. To circumvent the severe CFL constraint of explicit methods, implicit time dis- cretizations, which allow larger CFL numbers especially for stiff problems, are widely used in practice, especially in the CFD community to solve compressible flow prob- lems [27, 28, 6, 51, 50, 49, 46, 2]. Although most of the effort has been made for increasing accuracy of the time discretization and for increasing the efficiency of the nonlinear solver, there are not many works in the literature concerning the positivity- preserving property of implicit methods. For compressible turbulent flow problems, Batten et al. [3] have proposed a positive finite difference scheme by splitting the fluxes into “implicit” and “correction” parts and the source term into positive and negative parts. The “implicit” and negative terms are treated implicitly via the 15 Patankar trick [48]. In [43, 44, 45], Moryossef and Levy have developed implicit unconditional positive finite volume schemes for unsteady turbulent flows. Their main idea to preserve the positivity is to make the Jacobian matrix in each implicit time step an M -matrix. All these methods mentioned are low-order accurate and are complicated to generalize to high order. For DG methods, in [40], Meister and Ortleb have constructed an unconditional implicit positive DG scheme for solving shallow water equations. The positivity of the numerical approximation is preserved via a modified Patankar trick [48]. The method is shown to be conservative and uncondi- tional positivity-preserving, but only third-order accuracy is proved by a truncation error argument with no rigorous proof for arbitrary high-order spacial accuracy. In [82], Yuan, Cheng and Shu have developed a high-order unconditionally positive im- plicit DG method for radiative transfer equations. The positivity is preserved by utilizing the particular boundary conditions of the problem and by designing a novel rotational limiter. In this chapter, we extend the general framework in Chapter 1 to implicit tempo- ral discretizations and develop a positivity-preserving DG method with high-order spacial accuracy for one-dimensional conservation laws. We consider the DG method for the following two reasons. First, the discontinuous feature of its approximation space makes it a good fit for parallel implementation and for handling unstructured and extremely varying meshes. Second, for a class of implicit temporal discretiza- tions, including the backward Euler and Crank-Nicolson methods, Jiang and Shu have shown in [29] via the cell entropy inequality that the fully discrete schemes for the nonlinear conservation laws are unconditionally L2 -stable, which works for arbi- trary triangulation and any spacial order of accuracy. In this chapter, we consider the simplest implicit time discretization, i.e., the backward Euler temporal discretiza- tion as well as the Crank-Nicolson time discretization. Our focus is on constructing 16 a spatially high-order positivity-preserving DG scheme. The main conclusion is that in order to generalize the Zhang-Shu positivity-preserving limiter [86, 87] to the backward Euler DG scheme, a lower bound for the CFL number is required. This is proved theoretically for linear scalar equations and numerically verified for nonlinear equations. For the Crank-Nicolson method with P 1 element on uniform mesh, it is proved that both upper and lower bounds for the CFL number are required to have the positivity-preserving limiter take effect. The proposed positivity-preserving limiter is inexpensive and easy to implement. It not only preserves the positivity and high-order spacial accuracy but also makes the numerical scheme more robust, in the sense that it accelerates the convergence towards the steady-state solution and adds robustness to the nonlinear solver for extreme test cases. The organization of this chapter is as follows. In Section 3.2 the positivity- preserving technique is introduced. In particular, a CFL condition is derived to ensure the positiveness of the scheme. Then the introduction for the scheme for the compressible Euler system is given in Section 3.3. Numerical experiments are presented in Section 3.5 and concluding remarks are given in Section 3.6. 3.2 Positivity-preserving backward Euler DG scheme for scalar equations The backward Euler DG scheme for solving (3.1) is as introduced in Section 2.1, to seek un+1 h ∈ Vh such that in each subinterval Ij , we have (un+1 n+1 n h , v)j − ∆tLj (uh , v) = (uh , v)j (3.3) 17 for any v ∈ Vh . The spatial operator Lj is as defined in (2.5). The definition of a positivity-preserving DG scheme goes as below. Definition 3.2.1. A DG scheme is defined to be positivity preserving if given unh (x) ≥ 0, for any x ∈ Ω, then we have un+1 h (x) ≥ 0, ∀x ∈ Ω. In order to obtain more efficient implementations, this definition can be slightly relaxed to require positivity at specified quadrature points rather than at all the points. Definition 3.2.2 (Relaxed version). A DG scheme is defined to be positivity pre- serving if in each subinterval Ij , given unj (xαj ) ≥ 0, then we have un+1 j (xαj ) ≥ 0, where {xαj } is a set of selected quadrature points in cell Ij for each j = 1, · · · , N . In this section, we introduce how to add the positivity-preserving property to the scheme (3.3) by extending the framework in Section 2.2. As mentioned, the main difficulty lies in the first step, i.e., to derive sufficient conditions under which given unh (x) ≥ 0 for any x ∈ Ω, we have u¯n+1 j ∈ Ω for any j. For implicit DG methods, the polynomial un+1 j (x) depends on the previous information unh (x) in a global and implicit way. In this section, we would first show how to overcome this difficulty for scalar linear equations and derive a CFL condition. Then we will introduce the scaling limiter and then summarize the algorithm. 3.2.1 Preliminaries Let us first recall some definitions and results that will be useful in the following analysis. 18 Block circulant matrices For linear problems with periodic boundary conditions, when the mesh is uniform, the left hand side of the scheme (3.3) has a block circulant structure. The following block circulant matrices will become handy in the analysis for the uniform mesh cases. For a thorough introduction, one can refer to [15]. Definition 3.2.3. A matrix A, is defined to be a block circulant matrix bcirc(A1 , . . . , An ), if it has the following structure    A1 A2 ··· An    An A1 ··· An−1  A=      ···     A2 A3 ··· A1 The inverse of a block circulant matrix has the same structure. Lemma 3.2.1. The inverse of a block circulant matrix is block circulant. M-matrices Another useful tool is the so-called M -matrix. For a thorough introduction, one can refer to [4]. To define it let us first set Z n×n = {A = (aij ) ∈ Rn×n : aij ≤ 0, i 6= j} which is the set of all the n × n real matrices with nonpositive off-diagonal entries. In [52], the author listed forty equivalent characterizations for M -matrices. For our purpose, we adopt the following one as the definition. 19 Definition 3.2.4. A matrix A ∈ Z n×n is called an M -matrix if A is inverse-positive, that is, A−1 exists and each entry of A−1 is nonnegative. M -matrices have the following equivalent characterization [52]. Theorem 3.2.1. A matrix A = (aij ) ∈ Z n×n is an M -matrix if and only if aii > 0, 1 ≤ i ≤ n, and there exists a positive diagonal matrix D = diag{d1 , · · · , dn } such P that AD is strictly diagonally dominant, that is, aii di > j6=i |aij |dj for 1 ≤ i ≤ n. In particular, if D is the identity matrix, we have the following corollary. Corollary 3.2.1. A matrix A = (aij ) ∈ Z n×n is an M -matrix if aii > 0 and it is strictly diagonally dominant. Legendre polynomials In the following, we also utilize properties of Legendre polynomials. We consider the standard Legendre polynomials {pn (x)} defined on Iˆ by the following recursive relationship (n+1)pn+1 (x) = (2n+1)xpn (x)−npn−1 (x), p0 (x) = 1, p1 (x) = x, x ∈ Iˆ (3.4) In the following lemma, we collect some properties of Legendre polynomials that will be useful in the following analysis. For the proof, one can refer to [5]. Lemma 3.2.2. Legendre polynomials defined in (3.4) have the following properties 1. pn (1) = 1, pn (−1) = (−1)n . 2 R 2. p (x) pm (x) dx Iˆ n = δ , 2n+1 nm where δnm is the Kronecker delta. 20 d 3. (2n + 1)pn (x) = [p (x) dx n+1 − pn−1 (x)]. 1 dn 4. Rodrigues’ formula pn (x) = 2n n! dxn [(x2 − 1)n ]. 3.2.2 CFL condition for P 1 -DG methods on uniform meshes We consider the linear equation ut + ux = 0, x∈Ω (3.5) u(x, 0) = u0 (x), with periodic boundary condition. Then the scheme (3.3) becomes n+1 − − n+1 − + (un+1 n+1 n h , v)j − ∆t(uh , vx )j + ∆t[(uj+ 1 ) vj+ 1 − (uj− 1 ) vj− 1 ] = (uh , v)j (3.6) 2 2 2 2 Given unh (x) ≥ 0, ∀x ∈ Ω, we want to derive a CFL condition, under which the cell average u¯n+1 j ≥ 0, ∀j. For the linear equation (3.5) with linear DG methods (k = 1) on uniform meshes, the derivation of the CFL condition is straightforward due to the following reasons. 1. Given a linear function in Ij being positive is equivalent to say the function is positive at the two end points xj± 1 . And hence we can represent this polyno- 2 mial with the Lagrange interpolating polynomials at the two endpoints. 2. For uniform mesh, we are able to obtain an exact formula for the inverse operator and hence we can get a CFL condition that is necessary and sufficient for the positivity of the cell averages at the next time level. 21 If we choose a set of basis functions {φjl }kl=0 ∈ P k (Ij ), then (3.6) can be rewritten in the following matrix form M j cn+1 j − M j cnj + ∆tB1j cn+1 j − ∆tB2j cn+1 j n+1 j−1 − ∆tS cj =0 which is equivalent to (M j + ∆tB1j − ∆tS j )cn+1 j − ∆tB2j cn+1 j n j−1 = M cj (3.7) where Mlij = (φjl , φji )j , Slij = (∂x φjl , φji )j , (B1j )li = φjl (xj+ 1 )φji (xj+ 1 ), 2 2 (B2j )li = φjl (xj− 1 )φj−1 i (xj− 1 ). 2 2 If we transform all the integrations to the reference cell Iˆ via the affine mapping x−xj Tj (x) = hj /2 , we will obtain (M + 2λj B1 − 2λj S)cn+1 j − 2λj B2 cn+1 n j−1 = M cj (3.8) where Mli = (φˆl , φˆi )Iˆ, Sli = (∂x φˆl , φˆi )Iˆ, B1 = B1j , B2 = B2j The polynomials {φˆl }kl=0 are the basis of P k (I), ˆ with φlj = φˆl (Tj (x)). The coefficient vector cn+1 j is unchanged. For uniform mesh λj = λ, j = 1, · · · , N . If we set A = M + 2λB1 − 2λS and B = −2λB2 , the scheme (3.8) can be written in the following compact form Rcn+1 = Lcn =⇒ cn+1 = R−1 Lcn (3.9) 22 where R = bcirc(A, O, · · · , O, B) and L = Diag(M ). If we set  N ×(k+1)N ~ Θ1 0 · · · · · · 0  ~2 0 ··· 0    0 Θ Θ=      ···     ~ 0 · · · · · · 0 ΘN ~ j = (φ¯j0 , · · · , φ¯j ). where Θ k ¯ n+1 can be represented as in the following cell average Then the cell average u equation ¯ n+1 = ΘR−1 Lcn u (3.10) All we need to do is to work out the inverse matrix R−1 . In order to do this, first note that the matrix R a block circulant matrix and hence by Lemma 3.2.1, its inverse is still block circulant. We set R−1 = bcirc(D1 , D2 , · · · , DN ). Since RR−1 = I, direct calculation gives AD1 + BD2 = I, AD2 + BD3 = O AD3 + BD4 = O, ··· ADN −1 + BDN = O, ADN + BD1 = O And it is easy to obtain D1 = [A + (−1)N −1 B(A−1 B)N −1 ]−1 , (3.11) Dk = (−1)N +1−k (A−1 B)N +1−k D1 , k = 2, 3, · · · , N For the linear element (k = 1), we have the following result. Theorem 3.2.2. For backward Euler linear DG method, at time level n given 23 unh (x) ≥ 0 for any x ∈ Ω, then we have u¯n+1 j ≥ 0, for any j = 1, · · · , N if and only if 1 λ≥ (3.12) 3 Proof. We choose the basis {φˆ0 , φˆ1 } to be Lagrange interpolating polynomials at x = ±1. Then we have L = I and   2λ(3λ−1)  6λ2 +4λ+1 0 −A−1 B =   2λ(3λ+2) 6λ2 +4λ+1 0 2λ(3λ−1) 2λ(3λ+2) and if we set a = 6λ2 +4λ+1 and b = 6λ2 +4λ+1 , we furthermore have   N −1 a 0 (−A−1 B)N −1 =  aN −2 b 0 and C : = A + B(−A−1 B)N −1     N −1 λ + 1 −3λ   2λa 0 = +  λ 3λ + 1 −4λaN −1 0   N −1 λ + 1 + 2λa −3λ  =  λ − 4λaN −1 3λ + 1 It is easy to obtain det C = 6λ2 + 4λ + 1 − 2λ(3λ − 1)aN −1 = (6λ2 + 4λ + 1)(1 − aN ) Since |2λ(3λ − 1)| |a| = < 1, ∀λ > 0 6λ2 + 4λ + 1 24 we have det C > 0 and hence C is invertible and we have   1  3λ + 1 3λ D1 = C −1 =  det C −λ + 4λaN −1 λ + 1 + 2λaN −1   The local cell average operator θ~ := Θ ~ j = ( 1 , 1 ) for any j, and we claim that 2 2 ~ 1= 1 θD (2λ + 1 + 4λaN −1 , 4λ + 1 + 2λaN −1 )T > 0, ∀λ > 0 2 det C The second term is always positive since a ∈ (−1, 1). As for the first term, it is always positive when a > 0 or when N is an odd number. When N = 2m is an even number and a < 0, we have 8λ2 (3λ − 1) 2λ + 1 + 4λa2m−1 = 2λ + 1 + 4λ a a2m−2 ≥ 2λ + 1 + 4λa = 2λ + 1 + >0 6λ2 + 4λ + 1 for any λ > 0. Therefore, we have shown the previous statement. As for i > 1, we have Di = (−A−1 B)N +1−i D1   N +1−i a 0 =  D1 aN −i b 0   N +1−i N +1−i a (D1 )11 a (D1 )12  =  aN −i b(D1 )11 aN −i b(D1 )12 and then   ~ i= (D1 )11 N −i (D1 )12 N −i di := θD a (a + b), a (a + b) 2 2 25 Since 3λ + 1 3λ 12λ2 + 2λ (D1 )11 = > 0, (D1 )12 = > 0, a+b= > 0, ∀λ > 0 det C det C 6λ2 + 4λ + 1 we have 1 di > 0, ∀i ⇐⇒ a > 0 ⇐⇒ λ > 3 The cell average u¯n+1 1 is given by N X ~ˆ u¯n+1 1 = unh , φ) dj · (ˆ Iˆ j=1 Obviously, when λ > 31 , dj > 0 for any j and hence we have u¯n+1 j > 0. Otherwise, consider the following example, if we set uˆnh ≥ 0 to be supported only ~ˆ in the cell Ij where N − j is odd. Since moment vector (ˆ unh , φ) Iˆ is always positive, since for linear element the basis functions φˆl , l = 0, 1 are positive everywhere. When ~ˆ λ < 31 , the vector dj < 0, and hence u¯n+1 j = dj · (unh , φ) Iˆ < 0. Therefore, we have 1 u¯n+1 j ≥ 0, ∀j ⇐⇒ dj ≥ 0 ⇐⇒ λ ≥ 3 Remark 3.2.1. The sharpness of the bound can be verified by considering the fol- lowing example   1,  if x ∈ Ij u0h =  0, otherwise  After one step, we can check the minj u¯1j . 26 The analysis above, even though gives necessary and sufficient condition, relies on the uniform mesh assumption and works for small k. For non-uniform meshes, the matrix R is no longer a block circulant matrix and hence it is hard to obtain exact formula for R−1 . For large k, even for the uniform mesh case, the algebra will become very complicated and it seems impossible to extend to general polynomial degree k. 3.2.3 CFL condition for arbitrary polynomial degree on gen- eral meshes We continue considering the linear equation ut + ux = 0, x∈Ω (3.13) u(x, 0) = u0 (x), with periodic boundary condition. And still consider the following DG scheme n+1 − − n+1 − + (un+1 n+1 n h , v)j − ∆t(uh , vx )j + ∆t[(uj+ 1 ) vj+ 1 − (uj− 1 ) vj− 1 ] = (uh , v)j (3.14) 2 2 2 2 Now we consider the general mesh and arbitrary polynomial degree k. In this case, the inverse of the right hand side matrix can not be easily obtained. Therefore, we have to turn to other approaches. We first express u¯n+1 j in terms of unh . The idea is to take k + 1 different test functions as probes to extract the information out from un+1 n h (x) in terms of uh (x). First, let us take v = 1 in the scheme (3.14), we have u¯n+1 j + λj [(un+1 j+ 1 )− − (un+1 j− 1 )− ] = u¯nj , j = 1, · · · , N 2 2 27 ∆t where λj = hj . We can rewrite the system above in the matrix form as below Λ−1 u ¯ n+1 + A(un+1 )− = Λ−1 u ¯n (3.15) un1 , · · · , u¯nN )T , (un+1 )− = ((un+1 ¯ n = (¯ where u 3 )− , · · · , (un+1 N+ 1 )− ), Λ and A are N × N 2 2 matrices in the following form      λ1 0 ··· 0  1 0 · · · −1  ... .. .. .  . ..      0 λ2 .   −1 1 Λ= , A= .     .. .. .. .. . . ..   . . . 0     . . . 0      0 ··· 0 λN 0 · · · −1 1 Next, we need to express the cell boundary value (un+1 j+ 1 )− by unh . To this end, we 2 need to take other special test functions. Recall that the Dirac delta function has the following expansion ∞ 1X δ(x − y) = (2l + 1)pl (x)pl (y), x, y ∈ Iˆ 2 l=0 where pl (x) is the standard Legendre polynomial defined in (3.4). Then we set y = 1, truncate the summation at the (k + 1)th term and define k 1X k + 1 pk+1 (x) − pk (x) δˆk (x) = (2l + 1)pl (x) = , x ∈ Iˆ (3.16) 2 l=0 2 x−1 We have employed the Christoffel-Darboux formula [22] in the last equality. The following lemma says that this polynomial actually defines the Dirac delta ˆ at the point y = 1. function in P k (I) Lemma 3.2.3. The polynomial δˆk has the following properties 28 1. δˆk (x) ∈ P k (I) ˆ we have (w, δˆk ) ˆ = w(1). ˆ and for any w ∈ P k (I), I 2. In the cell Ij , define 2 ˆk δjk (x) = δ (Tj (x)) , x ∈ Ij (3.17) hj then it is the delta function in P k (Ij ) at xj+ 1 . 2 3. The mass is concentrated at x = 1, in the sense that for any k and j = 0, · · · , k − 1, we have (δˆk )(j) (1) − (δˆk )(j) (x) > 0 for any x ∈ [−1, 1). Proof. It is obvious that δˆk ∈ Pk (I). ˆ For any polynomial w ∈ Pk (I), ˆ we can write it as a linear combination of Legendre polynomials as below k X w(x) = cl pl (x), x ∈ Iˆ l=0 Then by the definition of δˆk and Lemma 3.2.2, we have k X k X k X (w, δˆk )Iˆ = cl (pl , δˆk )Iˆ = cl = cl pl (1) = w(1) l=0 l=0 l=0 The second property can be shown by a simple change of variable. For the third property, by the definition of δˆk , it suffices to show (i) (i) pl (1) − pl (x) > 0 (3.18) for l = 0, 1, · · · , k − 1, i = 0, . . . , l − 1 and for any x ∈ [−1, 1). It is straightforward to verify that (3.18) holds for l = 0, 1, 2. For l > 2, by Lemma 3.2.2, we have, (i) (i−1) (i) pl (x) = (2l − 1)pl−1 (x) + pl−2 (x) 29 If (3.18) holds for l ≤ m − 1, then we have when l = m, for fixed i < m and any x ∈ [−1, 1), (i−1) (i) (i−1) (i) p(i) (i) m (1) = (2m − 1)pm−1 (1) + pm−2 (1) > (2m − 1)pm−1 (x) + pm−2 (x) = pm (x) Then by induction, we have proved (3.18) and hence the third property in the lemma. With the help of the discrete delta function (3.17), we have the following repre- sentation of (un+1 j+ 1 )− . 2 Lemma 3.2.4. For linear scalar conservation law with f (u) = u discretized by the scheme (3.14), we have σjk (un+1 j+ 1 )− = ξjk u¯n+1 j unj , gjk )Iˆ + (ˆ (3.19) 2 where k−1 X k X σjk =1+ (2λj ) i+1 (αik − βik ), ξjk =2 (2λj )i βik , i=0 i=0 k−1 X h i gjk (x) = (2λj )i (δˆk )(i) (x) − βik i=0 and αik = (δˆk )(i) (1), βik = (δˆk )(i) (−1) (3.20) Proof. For fixed l = 0, 1, · · · , k − 1, take the test function to be (δjk )(l) (x) − (δjk )(l) (xj− 1 ) in the scheme (3.14). By the definition of δjk (x) and Lemma 3.2.3, 2 we have 30 (un+1 j , (δˆk )(l) (Tj (x)))j − hj βlk u¯n+1 j − 2λj (un+1 j , (δˆk )(l+1) (Tj (x)))j + ∆t(un+1 j+ 1 )− (αlk − βlk ) = (unj , (δˆk )(l) (Tj (x)))j − hj βlk u¯nj 2 If we expand un+1j in terms of the basis φlj (x) = pl (Tj (x)), l = 0, . . . , k, as un+1 j = Pk l n+1 l l=0 (cj ) φj and by a change of variable, we obtain un+1 (ˆ j , (δˆk )(l) )Iˆ − 2βlk u¯n+1 j un+1 − 2λj (ˆ j , (δˆk )(l+1) )Iˆ + 2λj (un+1 j+ 1 )− (αlk − βlk ) 2 unj , (δˆk )(l) )Iˆ − 2βlk u¯nj = (ˆ Pk where uˆnj = l n l=0 (cj ) pl (x). Or equivalently, un+1 (ˆ j , (δˆk )(l) )Iˆ − 2λj (ˆ un+1 j , (δˆk )(l+1) )Iˆ = 2βlk u¯n+1 j − 2λj (un+1 j+ 1 unj , (δˆk )(l) )Iˆ − 2βlk u¯nj )− (αlk − βlk ) + (ˆ 2 If we set un+1 Dl = (ˆ j , (δˆk )(l) )Iˆ, Cl = 2βlk u¯n+1 j − 2λj (un+1 j+ 1 unj , (δˆk )(l) )Iˆ − 2βlk u¯nj )− (αlk − βlk ) + (ˆ 2 then we have Dl − 2λj Dl+1 = Cl , l = 0, · · · , k − 1 and in particular, when l = k − 1, we have un+1 Dk−1 = 2λj Dk + Ck−1 = 2λj (ˆ j , (δˆk )(k) )Iˆ + Ck−1 = 4λj βkk u¯n+1 j + Ck−1 Then we have the following derivations (un+1 j+ 1 )− = (un+1 j un+1 , δjk )j = (ˆ j , δˆk )Iˆ = D0 = 2λj D1 + C0 2 31 l−1 X 2 l = (2λj ) D2 + 2λj C1 + C0 = · · · = (2λj ) Dl + Ci (2λj )i = · · · i=0 k−2 X = (2λj )k−1 [4λj βkk u¯n+1 j + Ck−1 ] + (2λj )i Ci i=0 k−1 X = 2(2λj )k βkk u¯n+1 j + (2λj )i Ci i=0 After plugging in the definition of Ci and after some manipulation, we obtain " k−1 # X 1+ (2λj )i+1 (αik − βik ) (un+1 j )− i=0 " k # k−1 X X = 2 (2λj )i βik u¯n+1 j + unj , (δˆk )(i) − βik )Iˆ (2λj )i (ˆ i=0 i=0 For the parameters {βjk , αjk }N j=1 , we have the following lemma, which will be useful in the proof of Proposition 3.2.1. Lemma 3.2.5. For any k ≥ 0, the following results hold 1. αik > βik for i = 0, · · · , k − 1 and hence σjk > 0 for any j. k k 2. βk−2i > 0, βk−2i−1 < 0, for i = 0, · · · , bk/2c. k k 3. βk−2i + βk−2i−1 > 0, for i = 0, · · · , bk/2c. Proof. The first statement is a direct conclusion of the third property of the discrete delta function in Lemma 3.2.3. For the second statement, let us first derive an explicit formula for βik . By the 32 Rodrigues’ formula in Lemma 3.2.2 we have 1 dl 2 l 1 dl  (x − 1)l (x + 1)l  pl (x) = l l (x − 1) = l l 2 l! dx 2 l! dx Then by the Leibniz’s rule, it is easy to obtain for i ≤ l,   (i) 1 i+l 1 (i + l)! l! pl (−1) = l l! [(x − 1)l ](i) |x=−1 = l (−2)l−i 2 l! l 2 i!l! (l − i)! (−2)−i (i + l)! = (−1)l i! (l − i)! and hence we have k 1X (−2)−i (l + i)! βik = (2l + 1) (−1)l 2 l=i i! (l − i)! k k 1 X (l + i)! l−i 1 X l = i+1 (2l + 1) (−1) = γ 2 i! l=i (l − i)! Ci l=i i (l+i)! where Ci = 2i+1 i! and γil = (2l + 1) (l−i)! (−1)l−i . For γil , we have (l + i + 1)! 2l + 3 l + i + 1 l |γil+1 | = (2l + 3) = |γ | > |γil |. (l + 1 − i)! 2l + 1 l − i + 1 i k Next for j = 0, · · · , bk/2c, let us consider βk−2j . If we replace i with k − 2j, we obtain k k k−1 k−2j−2 k−2j−1 k−2j Ck−2j βk−2j = (|γk−2j | − |γk−2j |) + · · · + (|γk−2j | − |γk−2j ) + |γk−2j | > 0. (3.21) k For βk−2j−1 , we have k k k−1 k−2j−2 k−2j−1 Ck−2j−1 βk−2j−1 = −[(|γk−2j−1 | − |γk−2j−1 |) + · · · + (|γk−2j−1 | − |γk−2j−1 )] < 0. (3.22) 33 k k Therefore, we can conclude that βk−2j > 0 and βk−2j−1 < 0. k k For the last statement, let us consider the sum βk−2j + βk−2j−1 . To this end, first for general γjl , let us consider the following expression 1 τjl := (|γjl | − |γjl−1 |) − |γj−1 l l−1 | + |γj−1 | 2j If we plug the definition of γjl in and after direct calculation, we obtain (l + j − 2)!  2 τjl = (l − j 2 )(2j + 1) + 2j − 1 ≥ 0.  j(l − j + 1)! If we combine (3.21) and (3.22) together we would obtain, j+1 X k−2j k k k−2i Ck−2j−1 (βk−2j + βk−2j−1 ) = τk−2j + |γk−2j |>0 i=0 which implies the desired conclusion. With Lemma 3.2.4 and the equation (3.15), we can obtain the following cell average equation. ¯ n+1 = L(un ) Tu (3.23) where   ξ1k 1 ξk σ1k + λ1 0 ··· − σNk  N   ξk ξ2k 1 .. ..   − σ1k σ2k + λ2 . .  T = 1   .. .. ..    . . . 0     ξk k ξN 1 0 ··· − σNk −1 k σN + λN N −1 and ! ! 1 gjk k gj−1 (L(un ))j = uˆnj , − k + uˆnj−1 , k , j = 1, . . . , N. 2λj σj σj−1 Iˆ Iˆ 34 A set of sufficient conditions to make the cell average u¯n+1 j positive for any j = 1, · · · , N are the following Condition I T is an M -matrix, or by Corollary 3.2.1 and Lemma 3.2.5, ξjk > 0, ∀j = 1, · · · , N (3.24) Condition II L(un ) is positive, or (ˆ unj−1 , gj−1 unj , 1/2λj − gjk /σjk )Iˆ + (ˆ k k /σj−1 )Iˆ > 0, j = 1, · · · , N (3.25) σjk By Lemma 3.2.3, it is easy to obtain gjk (x) < , ∀x ˆ for any k and j and hence ∈ I, 2λj the first term in (3.25) is always positive. Since uˆnj and uˆnj−1 can be any independent positive polynomials, the condition (3.25) is further reduced to (v, gjk )Iˆ ≥ 0, ˆ and v ≥ 0 ∀v ∈ P k (I) (3.26) i ˆk (i) Pk In summary, if we set Fk (λ, x) = i=0 (2λ) (δ ) (x) and then ξjk = Fk (λj , −1), gjk = Fk (λj , x) − Fk (λj , −1), sufficient conditions (3.24) and (3.25) actually require that the CFL number λj satisfy Fk (λj , −1) ≥ 0 (3.27) (Fk (λj , ·) − Fk (λj , −1), v(·))Iˆ ≥ 0, ˆ and v ≥ 0. ∀v ∈ P k (I) (3.28) The following theorem says that in order to make these two conditions hold simul- taneously, the CFL number λj can not be arbitrarily small. Theorem 3.2.3. When λj is small, conditions (3.27) and (3.28) can not hold at the same time. More specifically, we have the following two cases: 35 1. When the polynomial degree k is odd, there exists η1k > 0 such that when λj < η1k , the condition (3.27) does not hold. 2. When k ≥ 2 is even, there exists η2k > 0 such that when λj < η2k , the second condition (3.28) does not hold. Pk Proof. When k is odd, Fk (λj , −1) = i=0 βik (2λj )i is a polynomial of odd degree. By Lemma 3.2.5, we have the leading coefficient βkk > 0 and hence Fk (λj , −1) > 0 when λj is large. On the other hand, when λj = 0, Fk (0, −1) = β0k < 0, again by Lemma 3.2.5. Therefore, the polynomial Fk (λj , −1) must have at least one and at most k positive roots. If we take η1k to be the smallest one, we would have the first statement. When k ≥ 2 is even, in (3.28) take v = 1 and we obtain (Fk (λj , x) − Fk (λj , −1), 1)Iˆ k−1 X h i = i ˆ (2λj ) (δ )k (i−1) ˆk (i−1) (1) − (δ ) ˆk (i) (−1) − 2(δ ) (−1) ≥ 0 i=0 where (δˆk )(−1) (1) − (δˆk )(−1) (−1) := ˆk (x) dx. To simplify the notation, let us set R Iˆ δ y = 2λj and set k−1 X h i Gk (y) = i ˆk (i−1) y (δ ) ˆk (i−1) (1) − (δ ) ˆk (i) (−1) − 2(δ ) (−1) . i=0 Since k − 1 is odd, again, we would like to show it has positive real roots. First, when y = 0, we have 1−k Z k+1 1 Gk (0) = δˆk (x) dx − 2δˆk (−1) = 1 − = ≤− Iˆ 2 2 2 36 Next, let us check the leading coefficient of Gk (y). Since (δˆk )(k) = βkk , we have (δˆk )(k−2) = 21 βkk x2 + C1 x + C2 , where C1 and C2 are constants. As a consequence, we have k 1 1 αk−2 = βkk + C1 + C2 , k βk−2 = βkk − C1 + C2 , k βk−1 = −βkk + C1 2 2 and hence the leading coefficient of Gk (y) satisfies (δˆk )(k−2) (1) − (δˆk )(k−2) (−1) − 2(δˆk )(k−1) (−1) k k k = αk−2 − βk−2 − 2βk−1 = 2βkk > 0 Therefore, the odd-degree polynomial Gk (y) must have at least one and at most k −1 positive real roots. If we take η2k to be the smallest one, we can conclude the second statement. Theorem 3.2.3 indicates that unlike the situation for the DG method with Euler forward time discretization in [86], where an upper bound for the CFL number is sufficient for the cell average at the next time level to be positive, for the DG scheme with the Euler backward time discretization, the CFL number can not be arbitrarily small and an lower bound may be required. The following analysis confirms this statement. If we re-examine the condition (3.28) and note that for fixed λj , Fk (λj , x) ∈ ˆ the inner product in (3.28) can actually be approximated exactly by certain P k (I), Nq ˆ {ω α } the quadrature rules, say {(xα , ω α )}α=1 , where {xα } are the abscissas in I, weights and Nq is large enough such that the quadrature rule is exact for polynomials N of degree 2k. We denote {(xαj , ωjα )}α=1 q to be the transformed quadrature rule in Ij . 37 Then the condition (3.28) can be reduced to require Jαk (λj ) := Fk (λj , xα ) − Fk (λj , −1) ≥ 0, α = 1, · · · , Nq (3.29) If we define J0k (λj ) = Fk (λj , −1), together with condition (3.27), we require the CFL N number to make Nq +1 polynomials {Jαk (λj )}α=0 q positive. The following result states that it suffices to require λj ≥ 21 . Proposition 3.2.1. When minj λj ≥ 12 , we have Fk (λj , −1) > 0, (3.30) Fk (λj , x) − Fk (λj , −1) > 0, ∀x ∈ Iˆ (3.31) Proof. First, let us consider (3.30). By the definition we have Fk (λj , −1) = Pk k k i=0 (2λj ) βi , which can be rewritten as Fk (λj , −1) =  (2λj β2k + β1k ) + β0k , if k is even   (2λj )k−1 (2λj βkk + βk−1 k ) + ··· + (3.32) (2λj β1k + β0k ),   if k is odd When λj > 1/2, i.e., 2λj > 1, then by Lemma 3.2.5 it is easy to see that each term in (3.32) is strictly positive and hence Fk (λj , −1) > 0. For the condition (3.31), consider the ith derivative of Fk with respect to x. k ∂i (2λj )l−i (δˆk )(l) (x), X F k (λj , x) = i = 0, . . . , k ∂xi l=i 38 ∂k When i = k, by Lemma 3.2.3, we have F (λj , x) ∂xk k = βkk > 0 and hence ∂ k−1 ∂ k−1 k ˆ F k (λj , x) > Fk (λj , −1) = βk−1 + (2λj )βkk > 0, ∀x ∈ I. ∂xk−1 ∂xk−1 We have used the facts 2λj ≥ 1 and Lemma 3.2.5 to derive the last inequality. We can continue this procedure till i = 0 and obtain Fk (λj , x) > Fk (λj , −1), ∀x ∈ Iˆ Therefore, in summary, when minj λj ≥ 21 , we have conditions (3.30) and (3.31) hold and hence we can draw the conclusion. 1 Consequently, when λj ≥ 2 all the polynomials {Jαk }kα=0 are strictly positive. On the other hand by the proof of Theorem 3.2.3, at λj = 0 these polynomials can not be all nonnegative. This implies that if we denote Rk to be the set of all the positive roots of each polynomial Jαk , i.e., Rk = {r ∈ R+ : Jαk (r) = 0 for some α = 0, · · · , Nq } then we must have Rk 6= ∅, Rk ∈ (0, 21 ) and Rk is finite. Then if we set rk = max r, (3.33) r∈Rk we have the following theorem. Theorem 3.2.4. For the backward Euler DG scheme (3.14) with P k element for linear equation (3.13), at time level n given unj (xαj ) ≥ 0 at the quadrature points N {xαj }α=1 q ∈ Ij which give exact approximation for polynomials of degree 2k, we have 39 u¯n+1 j ≥ 0, for any j under the following CFL condition min λj ≥ rk (3.34) j where rk ∈ (0, 12 ) is defined as in (3.33). Proof. When λj ≥ rk by the definition of Rk , we have (3.29) and (3.27) hold. There- fore T in (3.23) is an M -matrix and L ≥ 0. By the definition of M -matrix, we can conclude the result. N Remark 3.2.2. We only require unh to be positive at quadrature points {xαj }α=1 q , which is weaker than the condition unh (x) ≥ 0, for any x ∈ Ω in the Definition 3.2.1. Remark 3.2.3. Even though the Theorem 3.2.4 is only proved for linear equations, numerical experiments suggest that for nonlinear equations, a lower bound for the CFL number is still necessary to make the cell average at the next time level positive. Remark 3.2.4. The results also hold for problems with positive source terms and positive inflow boundary conditions. Remark 3.2.5. The lower bound depends on the polynomial degree k as well as the quadrature rule we choose. For each k and fixed quadrature rule we can actually obtain the lower bound rk by solving the positive roots of each Jαk . In Table 3.1, we record the lower bounds for Legendre-Gauss-Lobbato (LGL) quadrature rule and the Legendre-Gauss (LG) rule respectively. We see that the LGL rule gives smaller lower bound. In practice we will limit the polynomial unj to make it positive at least at the N LGL points {xαj }α=1 q in each cell Ij . Remark 3.2.6. The lower bounds in Table 3.1 are sharp for odd k and sufficient 40 Table 3.1: Values of rk for k = 1, · · · , 5 and for Legendre-Gauss-Lobatto (LGL) and Legendre-Gauss quadrature rules respectively. LGL rule LG rule k Nq rk Nq rk 1 3 0.333 2 0.333 2 4 0.262 3 0.344 3 5 0.177 4 0.177 4 6 0.177 5 0.212 5 7 0.121 6 0.121 for even k. If we start from the following initial condition,   1,  if x ∈ IM u0 (x) = ,  0, otherwise  where M = arg maxj hj . After one step we record the minimum cell averages for different odd k and λ = rk ± 0.001 in Table 3.2. We see that when λ is slightly smaller than the lower bound rk , after one time step, at least one of the cell averages will become negative. If λ is larger than rk , the average will be uniformly positive. In Table 3.2: Minimum cell average after one time step for odd k. k λ = rk − 0.001 minj u¯1j λ = rk + 0.001 minj u¯1j 1 0.332 -3.498 E-04 0.334 5.284 E-35 3 0.176 -4.177 E-05 0.178 7.941 E-46 5 0.120 -2.135 E-06 0.122 1.980 E-59 1 the case k = 1, the lower bound r1 = 3 conincides with the lower bound in Theorem 3.2.2 for uniform meshes. 41 For even k, we consider a different initial condition   k  2(x−xM ) − 0.72 ,   hM if x ∈ IM u0 (x) =  0,  otherwise The Table 3.1 shows the minimum cell averages for different k and λ. We see that the lower bound for the CFL number is still necessary for the positivity of u¯1j and rk listed in Table 3.1 is sufficient. Table 3.3: Minimum cell average after one time step for even k. k λ minj u¯1j λ = rk minj u¯1j 2 0.170 -1.022 E-03 0.266 1.221 E-28 4 0.120 -1.343 E-03 0.177 5.021 E-46 3.2.4 Scaling limiter Once in each cell Ij , the cell average u¯n+1 j is positive, we limit the whole polynomial un+1 j (x) towards its cell average by utilizing the scaling limiter (2.12). In this case, we have G = [0, ∞) and ∂G = {0}. And then, the limiter can be defined in the following way, u˜n+1 j = θj [un+1 j − u¯n+1 j ] + u¯n+1 j , (3.35) where  ¯n+1 u if minx∈Ij un+1  j   n+1 ¯j −minx∈Ij u un+1 (x) , j (x) < 0 j θj = (3.36)  1,  otherwise Remark 3.2.7. In order to calculate the scaling parameter θj in (3.36) we need to calculate the minimum value of un+1 j in each cell Ij , which can be done efficiently up to k = 5 via the root formulas. For larger k, this calculation becomes expensive. But 42 recall that in Theorem 3.2.4, we only require each polynomial un+1 j to be positive on N LGL quadrature points {xαj }α=1 q and hence we can instead use the following scaling parameter  ¯n+1 u if minα un+1 (xαj ) < 0  j   ¯j −minα un+1 un+1 (xα , j j) θ˜j = j (3.37)  1,  otherwise in the limiter. This is similar with the scaling limiter in [86] and the same argument can be conducted here to prove that the modified limiter does not kill the original high-order accuracy either. 3.2.5 Algorithm for scalar equations Now, we can summarize the positivity-preserving algorithm for scalar equations as below. 1. At time level tn , given unj (x) being positive at least on the LGL quadrature N points {xαj }α=1 q . 2. Choose CFL number large enough a priori or adaptively enlarge it in each time step until u¯n+1 j ≥ 0, ∀j. 3. Apply the scaling limiter (3.35) with θj defined in (3.36) or (3.37) to un+1 j such that un+1 j (xαj ) ≥ , for any α = 1, . . . , Nq and j, where  is a small number to help get rid of the round-off effect. In the numerical examples, we take  = 10−13 . 43 3.3 Positivity-preserving backward Euler DG scheme for compressible Euler systems Next, let us consider the following compressible Euler system for ideal gas ut + f (u)x = 0, t ≥ 0, x ∈ [0, 1]. (3.38) with u = (ρ, m, E)T , f = (m, ρu2 + p, (E + p)u)T and 1 m = ρu, E = ρu2 + ρe, p = (γ − 1)ρe, 2 where ρ is the density, u the velocity, m the momentum, p the pressure, E the total energy and e is the internal energy. The constant γ > 1 is the ratio of specific heats. ∂f The Jacobian matrix ∂u has the following three eigenvalues ζ1 = u − c, ζ2 = u, ζ3 = u + c p where c = γp/ρ is the sound speed. The backward Euler DG scheme for (3.38) is defined to seek un+1 h ∈ Vh such that in each cell Ij we have (un+1 n+1 n h , v)j − ∆tLj (uh , v) = (uh , v)j , ∀v ∈ Vh (3.39) where Lj is defined as in (2.5) and the parameter in the global Lax-Friedrichs flux is set to be α = kc + |v|k∞ . 44 The physical solution lies in the following admissible set 1 m2     T G= u = (ρ, m, E) : ρ ≥ 0, p = (γ − 1) E − ≥0 (3.40) 2 ρ It can be verified that G is a convex set [87]. If we let unh to denote the DG approx- imation at time level n, then the positivity-preserving DG scheme for compressible Euler system is defined as below. Definition 3.3.1. A DG scheme for compressible Euler system (3.38) is positivity- preserving if at time level n given unh (x) ∈ G for all x ∈ Ω, then at the next time level (n + 1), we have un+1 h (x) ∈ G for all x ∈ Ω. Also, for the implementation reason, we have the relaxed version. Definition 3.3.2 (Relaxed version). A DG scheme for (3.38) is defined to be positivity preserving if in each subinterval Ij , given unj (xαj ) ∈ G, then we have un+1 j (xαj ) ∈ G, where {xαj } is a set of selected quadrature points in cell Ij for each j = 1, · · · , N . We design the positivity-preserving DG scheme by extending the general frame- work in Section 2.2. Unlike the scalar case, for system, we basically can not prove anything. But the analysis for the linear scalar equation suggests that starting from ¯ n+1 unj (x) ∈ G, in order to have the cell average u j ∈ G, a lower bound for the CFL number may be required. Therefore, we formulate the algorithm as follows. We stress on the applicability of the proposed algorithm rather than its theoretical justification. 1. At time level tn , in each cell Ij , given unj (xαj ) ∈ G at LGL quadrature points N {xαj }α=1 q . 45 2. Choose a large CFL number a priori or adaptively enlarge it in each time step ¯ n+1 until u j ∈ G for any j. 3. In each cell Ij , apply the following scaling limiter to the first component of N un+1 j to obtain ρn+1 j (x) ≥ 0 at LGL quadrature points {xαj }α=1 q ρ˜n+1 j = θ1 (ρn+1 j − ρ¯n+1 j ) + ρ¯n+1 j where  ρ¯n+1 , if minα ρn+1 (xαj ) < 0   j ρ¯j −minα ρn+1 n+1 j  (xα j) θ1 = j  1,  otherwise ˜ n+1 Denote the modified polynomial by u j . 4. In each cell Ij , apply the scaling limiter again to the whole modified polynomial ˜ n+1 u j such that p(xαj ) ≥ 0 for each α as below. n+1 n+1 n+1 u e un+1 e j = θ2 (e j −u ej ) + u ej where  p¯n+1 , if minα p˜n+1 (xαj ) < 0   j n+1 p¯j −minα p˜n+1 j  (xα j) θ2 = j  1,  otherwise and the pressure average is defined by p¯n+1 ¯˜ n+1 ). = p(u j j n+1 Remark 3.3.1. The modified polynomial u e ej ∈ G implemented in this way in general does not lies on ∂G. But this is a more robust version then the original one proposed in [87]. It was first developed in [72] and was further applied for the RHD system in [55]. The preservation of the high-order accuracy was verified via truncation error analysis [72] and numerical examples in [72, 55]. 46 3.4 Positivity-preserving Crank-Nicolson DG scheme The Crank-Nicolson P 1 -DG scheme for the conservation law (3.1) in each cell Ij , takes the following form ∆t (whn , v)j = (unh , v)j + Lj (unh , v), ∀v ∈ Vh 2 (3.41) ∆t (un+1 h , v)j − Lj (un+1 n h , v) = (wh , v)j , ∀v ∈ Vh 2 In this section we only consider P 1 element. We want to use the framework in Section 2.2 to construct positivity preserving DG scheme. The key is still to derive a CFL condition for the linear problem (3.13) under which given unh (x) ≥ 0, ∀x, we have u¯n+1 j ≥ 0, for any j. Note that in (3.41), the first step is a forward Euler scheme and the second step is a backward Euler scheme. A naive thinking is to first use the technique for the positivity-preserving forward Euler scheme to derive a CFL condition, under which whn ≥ 0 at least at selected quadrature points. Then once whn (x) ≥ 0, we can use the result for the backward Euler scheme in the last section to obtain u¯n+1 j ≥ 0. Unfortunately, the following analysis shows the CFL condition for the forward Euler step does not exist even for P 1 element. ˆ Set the discrete delta Let us consider any xαj ∈ Ij and set xα = Tj (xαj ) ∈ I. ˆ as δˆα1 , which by function in P 1 (Ij ) at xαj as δα1 and the corresponding one in P 1 (I) definition takes the following form 1 3xα δˆα1 = + x 2 2 If we take v = δα1 in the forward Euler scheme in (3.41) and after a change of variable, 47 we obtain unh , δˆα1 ))Iˆ − 2λj [(unj+ 1 )− δˆα1 (1) − (unj− 1 )− δˆα1 (−1)] wαn = unα + 2λj (ˆ 2 2 where wαn = wˆjn (xα ) = wjn (xαj ) and the same for unα . After plugging the definition of δˆα1 into the scheme, we obtain wαn = unα + 6λj xα u¯nj − λj unh (x− j+ 1 )(1 + 3xα ) + λj unh (x− j− 1 )(1 − 3xα ) 2 2 = [unα − λj unj (xj+ 1 ) + 3λj xα unj (xj− 1 )] + λj (1 − 3xα )unj−1 (xj− 1 ) 2 2 2 Since unj and unj−1 are independent positive polynomials, if xα > 13 , the we can set unj−1 (xj− 1 ) arbitrarily large and hence there exists no λj making wαn positive. Recall 2 in Table 3.1, for k = 1, we need to require wjn to be positive at at least three LGL points and two GL points, which all include abscissas greater than 13 . Therefore, we can not combine the results for the positivity-preserving forward and backward Euler methods and extend it to the Crank-Nicolson method. 3.4.1 CFL condition for P 1 -DG on uniform meshes If we further assume the mesh is uniform, by brutal force as for the backward Euler method, we can get some necessary and sufficient results for the Crank-Nicolson method. Similar with (3.7), the Crank-Nicolson time discretization for the linear DG method for the linear problem (3.13), which, in each cell Ij can be written in the 48 following matrix form. cn+1 j − cnj 1  j n+1 j (S cj − B1j cn+1 + B2j cn+1 j n j n j n  M = j j−1 ) + (S cj − B1 cj + B2 cj−1 ) ∆t 2 which is equivalent to   ∆t j ∆t j n+1 ∆t j n+1 ∆t j n j M − S + B1 cj − B2 cj−1 = Mj cnj + (S cj − B1j cnj + B2j cnj−1 ) 2 2 2 2 After rescaling, we have (M − λj S + λj B1 )cn+1 j − λj B2 cn+1 n n j−1 = (M + λj S − λj B1 )cj + λj B2 cj−1 (3.42) or in the compact form P cn+1 = P˜ cn where for uniform mesh with λ := λ1 = · · · = λN , we denote P = bcirc{M − λS + λB1 , O, · · · , O, −λB2 } and P˜ = bcirc{M + λS − λB1 , O, · · · , O, λB2 }. Then we have the following cell average equation ¯ n+1 = ΘP −1 P˜ cn u (3.43) For the linear element, we have the following result. Theorem 3.4.1. For the Crank-Nicolson DG method with P 1 element for the linear problem (3.13), given unh (x) > 0 for any x ∈ Ω, we have u¯n+1 j > 0 for any j, if and only if   2 λ∈ , λN (3.44) 3 q 2 with λN & 3 as N → ∞. 49 λ Proof. Note that if we replace λ with 2 in the scheme (3.8) the inverse of P takes exactly the same form as R−1 , i.e., P −1 = bcirc{D1 , D2 , · · · , DN }, with Di defined in (3.11). If we set A˜ = M + λ(S − B1 ), ˜ = λB2 B then we have P −1 P˜ = bcirc{D1 A˜ + D2 B, ˜ D2 A˜ + D3 B, ˜ · · · , DN A˜ + D1 B}. ˜ When i = 2, · · · , N − 1, we have αλ(9λ2 − 4)   ~ i A˜ + Di+1 B) di = θ(D ˜ = N −i , 3λαa (a + b) 3λ2 + 4λ + 2 where α = det C is always positive. Then 2 di ≥ 0, ∀i = 2, · · · , N − 1 ⇐⇒ 9λ2 − 4 > 0 & aN −i > 0, ∀i ⇐⇒ λ > 3 When i = N , similar procedure gives   6λ(3λ + 1)(λ + 4) dN = , 9αλ(a + b) > 0, ∀λ > 0. 3λ2 + 4λ + 2 At last, let us examine the case when i = 1. In this case, direct calculation gives d1 = (pN (λ), qN (λ)) 50 where 3 pN (λ) = − λ2 (1 − aN −2 b) + λaN −2 (a − b) + 1 2 3 qN (λ) = − λ2 (1 − aN −1 ) + λ(2 + aN −1 ) + 1 2 Let us first examine the first entry pN (λ). We claim it has the following three properties. (P1) pN (λ) is strictly decreasing and has unique positive root λN > 0. (P2) For fixed λ > 23 , we have pN +1 (λ) < pN (λ). (P3) For each N , λN > 32 . The proof for these statements are postponed after the proof of this theorem. With q these statements in hand, we want to show λN & 23 . First of all, since |a| < 1, for fixed λ > 0, it is easy to see that as N → ∞ 3 pN (λ) → p0 (λ) = − λ2 + 1 2 q 2 which has positive root λ0 = 3 . By (P2), we have for any M > N , pN (λ0 ) > pM (λ0 ). Then by letting M → ∞, we will have pN (λ0 ) > p0 (λ0 ) = 0. Then by (P1), we have λN > λ0 . Therefore, {λN } is a decreasing sequence with lower bound λ0 and ˜ 0 = limN →∞ λN . Next, we want to show λ hence has limit λ ˜ 0 = λ0 . Suppose λ ˜ 0 > λ0 , ˜ 0 ) < 0. Since pN (λ then p0 (λ ˜ 0 ) → p0 (λ ˜ 0 ) there must exist N0 such that pN (λ ˜0) < 0 0 q and hence λN0 < λ ˜ 0 which is a contradiction. Therefore, λN & λ0 = 2 . 3 Then let us look at qN (λ). Similar with pN , one can show the following results. (Q1) When λ > 23 , qN (λ) > qN +1 (λ). 51 2 0 (Q2) When λ > 3 and N ≥ 3, we have qN (λ) < 0. (Q3) When λ > 23 , qN has unique positive root rN . √ Then by the same token, one can show that rN & r0 = (2 + 10)/3. In summary, for the Crank-Nicolson DG with linear element, the scheme is mono- tone if and only if   2 λ∈ , λN 3 q 2 with λN & 3 . Remark 3.4.1. The sharpness of the bounds can be verified by considering the fol- lowing two examples. The first one is for testing the upper bound. Consider the following initial condition.   (x − x  j− 12 )/h, if x ∈ Ij u0h = (3.45)  0,  otherwise After one step, we can check the minj u¯1j . In Table 3.4, we list the minimum cell q q average for λ = 3 −10 and λ = 23 +10−9 respectively. We see that for small N , 2 −9 q when λN > 23 , we have the positive minimum cell averages after one time step. As p we increase the number of cells N , the upper bound converges to 2/3 and hence the minimum cell average becomes negative. This verifies the dependence of the upper q 2 bound on the cell number N . But if we take λ slightly smaller than 3 the cell average will always be positive. 52 Table 3.4: Values of minj u¯1j with initial condition (3.45) for N = 4, 8, 16, 32 and different λ. q minj u¯1j q N λ = 23 + 10−9 λ = 23 − 10−9 4 6.12 e-5 6.12 e-5 8 6.11 e-11 7.35 e-10 16 -3.37 e-10 3.34 e-19 32 -3.37 e-10 5.99 e-40 For the lower bound, let us consider the following initial condition.   (−x + x  j+ 21 )/h, if x ∈ Ij u0h = (3.46)  0,  otherwise In Table 3.4, we list the same results for λ = 32 ± 10−9 . We see that when λ is slightly smaller than the lower bound, the minimal cell average is negative. However, when λ is slightly larger than 32 , as expected, the cell average is always positive. And the lower bound is independent of the cell number N . Table 3.5: Values of minj u¯1j with initial condition (3.46) for N = 4, 8, 16, 32 and different λ. minj u¯1j N λ = − 10−9 λ = 23 + 10−9 2 3 4 -1.48 e-10 8.53 e-18 8 -1.48 e-10 9.51 e-28 16 -1.48 e-10 2.54 e-119 32 -1.48 e-10 9.51 e-28 Remark 3.4.2. The CFL condition derived in 3.4.1 is very stringent, which makes tthe positivity-preserving Crank-Nicolsono DG methods impractical. In the numerical experiments, we only implement the Crank-Nicolson P 1 DG method for accuracy test purpose, without any other applications. 53 Remark 3.4.3. The conclusion for Crank-Nicolson is surprising at the first look. But if we recall the relationship between C-N and forward and backward Euler meth- ods, this actually is a very natural result. Here is the heuristic argument. The Crank-Nicolson method for ut = L(u) can be written as u(1) = un + ∆tL(un )   n+1 1 (1) 1 n ∆t u = u + u + L(un+1 ) 2 2 2 which is actually a convex combination of an Euler forward and an Euler backward methods. Then u(1) ≥ 0 requires λ ≤ rN λ 1 2 un+1 ≥ 0 requires > ⇐⇒ λ > 2 3 3 and hence we have to require both a lower and a upper bound for the Crank-Nicolson to be positivity-preserving. Proof of (P1) Let us prove that pN is a strictly decreasing function for each N . Consider the derivative of it. Direct calculation gives 3λα2 p0N (λ) = gN (λ) + 2, βN where α = 3λ2 − 2λ, β = 3λ2 + 4λ, gN (λ) = αN −2 9λ4 + (9N + 6)λ3 + (24N − 20)λ2 + (10N − 16)λ − 4N ] − β N .  54 We want to show gN (λ) < 0 for any λ > 0 and N . First of all, it is very easy to check the result for N = 1, 2, 3 and hence we assume N ≥ 4. And discuss in the following several situations. First, when αN −2 < 0, i.e., when λ < 2 3 and N is odd, since in this case, |α| < 1, we have gN ≤ 4N |α|N −2 − β N ≤ 4N − β N ≤ 4N − 2N − 3N λ2N < 0, (3.47) for N ≥ 4. Next, when αN −2 ∈ (0, 1), in this case we must have λ < 1 and hence gN ≤ 9λ4 + (9N + 6)λ3 + (24N − 20)λ2 + (10N − 16)λ − β N ≤ 9λ4 + (9N + 6)λ3 + (24N − 20)λ2 − (4λ + 2)N N X 4 3 2 = 9λ + (9N + 6)λ + (24N − 20)λ + (10N − 16)λ − Cni (4λ)i 2N −i i=0 N X =: ci λ i i=0 It is obvious that ci < 0 for i = 0 and i > 4. Let us check for ci when i = 1, 2, 3, 4 and for N ≥ 4 1. When i = 1, c1 = 10N − 16 − N − N 2N +1 < 0. 2. When i = 2, we have c2 = 24N − 20 − CN2 42 2N −2 = 24N − 20 − 2N +3 (N 2 − N ) ≤ 24N − 20 − 23 (N 2 − N ) = −8N 2 + 32N − 20 < 0, 55 3. When i = 3, we have c3 = 9N + 6 − CN3 2N +3 = 9N + 6 − 2N +3 (N 3 − 3N 2 + 2N ) = −2N +3 (N 3 − 3N 2 ) − (2N +4 N − 9N − 6) <0 4. For c4 , we have c4 = 9 − CN4 44 22 < 0. Therefore, gN < 0 when αN −2 ∈ (0, 1). At last, let us consider the case where αN −2 > 1, i.e., when λ > 1. First, we have gN ≤ (3λ2 )N −2 [9λ4 + (9N + 6)λ3 + (24N − 20)λ2 + (10N − 16)λ] − (3λ2 + 4λ)N ≤ 3N λ2N + 3N −1 (3N + 2)λ2N −1 + 3N −2 (24N − 20)λ2N −2 + 3N −2 (10N − 16)λ2N −3 N X − Cnk 3N −i 4i λ2N −i i=0 XN =: di λi i=0 Again, we check dN −i for i = 0, 1, 2, 3 and for N ≥ 4. 1. Consider i = 0, dN = 3N − 3N = 0. 2. For i = 1, we have dN −1 = 3N −1 (3N + 2) − 4N 3N −1 = 3N −1 (3N + 2 − 4N ) = −3N −1 (N − 2) < 0. 3. For i = 2, we have dN −2 = 3N −2 (−8N 2 + 32N − 20) < 0. 56 4. For i = 3, we have   N −3 32N (N − 1)(N − 2) dN −3 = 3 3(10N − 16) − 3 = 3N −4 [90N − 144 − 32N (N − 1)(N − 2)] ≤ 3N −4 N (−32N 2 + 96N + 26) <0 In conclusion, we have shown gN < 0, which immediately implies p0N < 0, for any λ > 0. Proof of (P2) For fixed λ > 32 , we have 1 3 (pN +1 (λ) − pN (λ)) = b λ2 (a − 1) + λ(b − 1)(a − 1) aN −2 2 3bλ2   = (a − 1) + (a − b)λ 2 3λ2 (3λ − 2)(λ + 2) = (a − 1) 2β <0 since |a| < 1. Proof of (P3) Because of (P1), it suffices to show pN (2/3) > 0 and pN (2) < 0. This is straightforward as below.  2 3 2 pN (2/3) = − >0 2 3 1  N −2 48 4 48 pN (2) = −5 + ≤ −5 + < 0, ∀N ≥ 2. 11 11 11 57 3.5 Numerical experiments In this section, we present numerical examples. First, we verify the accuracy of the Crank-Nicolson positivity-preserving P 1 -DG method on linear problem. Then for the backward Euler DG method, we verify the high-order spacial accuracy of the proposed method by testing it on both linear and nonlinear steady-state problems. An acceleration of the convergence towards the steady state solution is also observed. Next, we test the methods on moving-shock problems. At last, examples for com- pressible Euler system will be presented. In all the examples, the domain is first uniformly decomposed and then each node is randomly perturbed up to 20% of the mesh size. To further solve the nonlinear system (3.3) and (3.39), there have been many works on how to build efficient solvers, such as the work in [51, 50, 49]. But since our main focus is on the positivity preserving property rather than the efficiency of the nonlinear solver, we use the Newton method [16] for the nonlinear system up to accuracy 10−13 . For the robustness and accuracy reasons, in each Newton iteration, the Jacobian matrix is solved with the direct solver. 3.5.1 Accuracy tests First, let us test the accuracy of the proposed method. The accuracy is first tested on the linear problem solved with the Crank-Nicolson P 1 -DG method. And then we check the high-order spacial accuracy with the steady-state solution to both of the linear equation and the Burgers equation. For the steady-state simulation, we take ∆t = 10 maxj hj and march in time until kun+1 h − unh k2 ≤ 10−12 . 58 Example 3.5.1.1 (Linear problem with Crank-Nicolson P 1 -DG). Consider the fol- lowing problem. ut + ux = 0, x ∈ [0, 2π] u0 (x) = sin4 (x) with periodic boundary condition. The exact solution is u(x, t) = sin4 (x − t). The equation is solved with P 1 -DG with Crank-Nicolson time discretization. Uniform mesh is employed. As the Theorem 3.4.1 indicates, the positivity-preserving limiter q will only work when λ ∈ [ 32 , 23 ]. Therefore, we take λ = 0.81 in this example. In Table 3.6, we record numerical results at T = 1. The L2 and L∞ errors and the numerical orders of convergence as well as the minimal values at time T are presented. We see that without the positivity-preserving limiter, the minimum value are negative. With the limiter on, the positivity of the numerical solution is preserved and the numerical order of accuracy is maintained. Table 3.6: Error table for Example 3.5.1.1 at T = 1 with λ = 0.81, uniform mesh. N L2 error L∞ error min uh 40 4.91 E-2 – 9.11 E-2 – -3.48 E-2 80 1.31 E-2 1.91 2.42 E-2 1.91 -7.93 E-3 no limiter 160 3.34 E-3 1.97 6.08 E-3 1.99 -1.37 E-3 320 8.38 E-4 1.99 1.52 E-3 2.00 -2.29 E-4 40 4.74 E-2 – 9.11 E-2 – 4.96 E-4 80 1.34 E-2 1.82 2.63 E-2 1.79 2.10 E-5 with limiter 160 3.39 E-3 1.99 6.08 E-3 2.11 1.55 E-6 320 8.42 E-4 2.01 1.52 E-3 2.00 1.00 E-7 Example 3.5.1.2 (Steady-state solution to linear problem). For the linear equation, we consider the following problem ut + ux = sin4 (x) u(x, 0) = sin2 (x) (3.48) u(0, t) = 0 59 with the outflow boundary condition at x = 2π. The exact solution u(x, t) can be derived by the characteristic theory and can be shown to be positive for all t > 0. In Table 3.7 and Table 3.8 we record the errors, numerical orders of accuracy and the minimum value of the numerical approximation, with and without the positivity- preserving limiter respectively. We see that without the positivity-preserving limiter the minimum value of the steady-state approximation is negative. When the limiter is put on, the minimum value becomes positive and the high-order accuracy is not destroyed. Table 3.7: Error table for Example 3.5.1.2, approximation of the steady state solution to the linear problem (3.48), without the positivity preserving limiter. k N L2 error order L∞ error order min uh 20 4.555 E-2 – 6.376 E-2 – -6.037 E-2 40 1.177 E-2 1.95 1.704 E-2 1.90 -5.075 E-3 1 80 2.967 E-3 1.99 4.347 E-3 1.97 -2.808 E-4 160 7.434 E-4 2.00 1.092 E-3 1.99 -1.170 E-5 320 1.859 E-4 2.00 2.733 E-4 2.00 -3.901 E-7 20 5.749 E-3 – 9.561 E-3 – -1.667 E-3 40 7.482 E-4 2.94 1.121 E-3 3.09 -6.550 E-5 2 80 9.449 E-5 2.99 1.543 E-4 2.86 -2.163 E-6 160 1.184 E-5 3.00 1.975 E-5 2.97 -6.853 E-8 320 1.481 E-6 3.00 2.484 E-6 2.99 -2.149 E-9 20 6.987 E-4 – 6.013 E-4 – -8.652 E-4 40 4.564 E-5 3.94 4.743 E-5 3.66 -3.909 E-5 3 80 2.885 E-6 3.98 2.986 E-6 3.99 -1.329 E-6 160 1.808 E-7 4.00 1.887 E-7 3.98 -4.240 E-8 320 1.131 E-8 4.00 1.178 E-8 4.00 -1.332 E-9 20 6.094 E-5 – 4.622 E-5 – -6.545 E-6 40 1.955 E-5 4.96 1.335 E-6 5.11 -1.329 E-6 4 80 6.149 E-8 4.99 4.597 E-8 4.86 -5.211 E-8 160 1.925 E-9 5.00 1.471 E-9 4.97 -1.715 E-9 320 6.018 E-11 5.00 4.623 E-10 4.99 -5.426 E-11 Example 3.5.1.3 (Steady-state solution to Burgers problems). For the Burgers 60 Table 3.8: Error table for Example 3.5.1.2, approximation of the steady state solution to the linear problem (3.48) with the positivity preserving limiter. k N L2 error order L∞ error order min uh 20 4.253 E-2 – 6.376 E-2 – 1.000 E-13 40 1.173 E-2 1.86 1.704 E-2 1.90 1.000 E-13 1 80 2.966 E-3 1.98 4.347 E-3 1.97 1.000 E-13 160 7.434 E-4 2.00 1.092 E-3 1.99 1.000 E-13 320 1.859 E-4 2.00 2.733 E-4 2.00 1.000 E-13 20 5.762 E-3 – 9.561 E-3 – 1.000 E-13 40 7.482 E-4 2.95 1.121 E-3 3.09 1.000 E-13 2 80 9.449 E-5 2.99 1.543 E-4 2.86 1.000 E-13 160 1.184 E-5 3.00 1.975 E-5 2.97 1.000 E-13 320 1.481 E-6 3.00 2.484 E-6 2.99 1.000 E-13 20 1.015 E-3 – 2.240 E-3 – 1.000 E-13 40 5.077 E-5 4.32 9.673 E-5 4.53 1.000 E-13 3 80 2.932 E-6 4.11 3.253 E-6 4.89 1.000 E-13 160 1.812 E-7 4.02 1.887 E-7 4.11 1.000 E-13 320 1.131 E-8 4.00 1.178 E-8 4.00 1.000 E-13 20 6.141 E-5 – 4.622 E-5 – 1.000 E-13 40 2.230 E-5 4.78 4.356 E-6 3.41 1.000 E-13 4 80 6.830 E-8 5.03 1.700 E-7 4.68 1.000 E-13 160 2.045 E-9 5.06 5.591 E-9 4.93 1.000 E-13 320 6.214 E-10 5.04 1.773 E-10 4.98 1.000 E-13 equation, we consider the steady state solution to the following problem u2   x ut + = sin 2 x 4 u(x, 0) = x (3.49) u(0, t) = 0 with the outflow boundary condition at x = 2π. Again, by the characteristic theory, one can show that the solution to this problem is always positive. In Table 3.9 and 3.10, we present numerical results for the cases where the positivity-preserving limiter is on and off respectively. This example shows the effectiveness of the positivity preserving limiter for nonlinear scalar problem. With the limiter on, the solution 61 stays positive and the high-order accuracy is preserved. Table 3.9: Error table for Example 3.5.1.3, approximation of the steady state solution to the Burgers equation (3.49), without the positivity preserving limiter. k N L2 error order L∞ error order min uh 20 2.915 E-6 – 5.085 E-6 – -8.073 E-6 40 3.508 E-7 3.05 6.358 E-7 3.00 -1.009 E-6 2 80 4.293 E-8 3.03 7.948 E-8 3.00 -1.261 E-7 160 5.304 E-9 3.02 9.934 E-9 3.00 -1.577 E-8 20 2.336 E-9 – 1.648 E-9 – -3.855 E-10 40 1.445 E-10 4.02 1.030 E-10 4.00 -1.204 E-11 3 80 9.010 E-11 4.00 6.525 E-11 3.98 -3.763 E-13 160 6.031 E-12 3.90 4.978 E-12 3.71 -1.175 E-14 20 3.403 E-10 – 6.312 E-10 – -1.497 E-9 4 40 1.010 E-11 5.07 1.970 E-11 5.00 -4.678 E-11 80 3.116 E-13 5.02 5.398 E-13 5.19 -1.462 E-12 Table 3.10: Error table for Example 3.5.1.3, approximation of the steady state solu- tion to the Burgers equation (3.49), with the positivity preserving limiter. k N L2 error order L∞ error order min uh 20 3.207 E-6 – 4.408 E-6 – 1.000 E-13 40 3.696 E-7 3.12 6.010 E-7 2.87 1.000 E-13 2 80 4.435 E-8 3.06 8.171 E-8 2.88 1.000 E-13 160 5.414 E-9 3.03 1.083 E-8 2.92 1.000 E-13 20 2.338 E-9 – 1.648 E-9 – 1.000 E-13 40 1.445 E-10 4.02 1.030 E-10 4.00 1.000 E-13 3 80 9.010 E-11 4.00 6.511 E-11 3.98 1.000 E-13 160 6.031 E-12 3.90 4.987 E-12 3.71 1.051 E-14 20 4.792 E-10 – 9.061 E-10 – 1.000 E-13 4 40 1.253 E-11 5.26 3.102 E-11 4.87 1.000 E-13 80 3.653 E-13 5.10 1.086 E-12 4.84 1.000 E-13 62 Next, let us turn to another steady-state Burgers problem u2   x ut + = sin3 2 x 4 x u(x, 0) = sin2 (3.50) 4 u(0, t) = 0 with the outflow boundary condition at x = 2π. This problem is more difficult than the previous one, since the source is much closer to zero around x = 0. In Table 3.11, we present the error, numerical convergence rate as well as the number of time steps taken to reach the steady state solution, which is denoted by NT . Both cases where the positivity preserving limiter is on and off are presented respectively. We use this example to illustrate that the stability added to the scheme by the positivity- preserving limiter helps accelerate the convergence towards the steady-state solution. Table 3.11: Error table for Example 3.5.1.3, approximation of the steady-state so- lution to the Burgers equation (3.50), with and without the positivity preserving limiter. without limiter with limiter 2 2 k N L error order NT L error order NT 20 1.71 E-5 – 670 1.71 E-5 – 158 40 2.11 E-6 3.02 1962 2.11 E-6 3.02 500 2 80 2.62 E-7 3.01 4973 2.62 E-7 3.01 1560 160 3.29 E-8 2.99 8352 3.27 E-8 3.01 4560 20 9.10 E-8 – 1162 9.10 E-8 – 970 3 40 5.36 E-9 4.08 3440 5.36 E-9 4.08 2211 80 3.44 E-10 3.96 9007 3.35 E-10 4.00 2653 63 3.5.2 Moving shocks Next, we test the proposed scheme on problems involving moving shocks. In all the numerical experiments below, quadratic P 2 element and non-uniform mesh are employed. Example 3.5.2.1 (Linear Problem). The first example is the linear equation with the initial data   1,  if x ∈ [3, 4] u0 (x) =  0, otherwise  In Figure 3.1, we present the numerical results at T = 0.2, with and without the positivity-preserving limiter respectively. For the linear equation, we take minj λj = r2 = 0.266 as listed in Table 3.1. As the zoom-in plots show, the limiter helps the solution keep positive. Without the limiter a negative undershoot appears. Example 3.5.2.2 (Burgers problem). Next, let us consider the Burgers equation with the initial condition u0 (x) = 1 + sin(x) and periodic boundary conditions. The initial condition is positive and takes value zero at x = π. At T = 1.5, a shock is developed. In Figure 3.2, we show the numerical approximation with and without the limiter respectively. We see that even though the profiles are smeared due to the first-order accuracy of the backward Euler time discretization, the effectiveness of the positivity-preserving limiter can still be observed in the zoom-in plots around the shock. Example 3.5.2.3 (Buckley-Leverett problem). In this example, the Buckley- 4u2 Leverett problem is considered, with f (u) = 4u2 +(1−u)2 . For this flux function, we have f 0 (0) = f 0 (1) = 0. If we start with the usual step function u0 (x) = I[3,4] , numer- ical experiments indicate that no matter how large ∆t is, we can not make the cell average u¯1j positive for all j. The same phenomena is also observed for the Burgers 64 0.1 1 0.8 0.05 0.6 u u 0 0.4 0.2 -0.05 0 -0.1 1 2 3 4 5 6 1 1.5 2 2.5 3 x x (a) without limiter (b) zoom in around x = 2 0.1 1 0.8 0.05 0.6 u u 0 0.4 0.2 -0.05 0 -0.1 1 2 3 4 5 6 1 1.5 2 2.5 3 x x (c) with limiter (d) zoom in around x = 2 Figure 3.1: Example 3.5.2.1: at T = 0.2, with ∆t = 0.266 maxj hj and N = 120. Solid line is the exact solution. Blue squares are point values on the Legendre-Gauss- Lobatto points. problem with step function as the initial condition. This indicates that the lower bound for the CFL number is also necessary for the positivity-preserving limiter to work for nonlinear problems. And in general, the limiter may not be applied for problems with sonic points. In this example, we change the initial condition to   0.9,  if x ∈ [3, 4] u0 (x) = 10−3 , otherwise   65 0.3 2 0.25 0.2 1.5 0.15 u 1 u 0.1 0.05 0.5 0 -0.05 0 -0.1 0 1 2 3 4 5 6 4 4.2 4.4 4.6 4.8 5 5.2 5.4 x x (a) without limiter (b) zoom in around the shock 0.3 2 0.25 0.2 1.5 0.15 u 1 u 0.1 0.05 0.5 0 -0.05 0 -0.1 0 1 2 3 4 5 6 4 4.2 4.4 4.6 4.8 5 5.2 5.4 x x (c) with limiter (d) zoom in around the shock Figure 3.2: Example 3.5.2.2 at T = 1.5. Solid line is the exact solution. Blue squares are point values of the numerical approximation on the Legendre-Gauss- Lobatto points. Mesh with N = 120 cells and CFL = 2. In Figure 3.3, we present the plots with and without limiter. Even though the profile is smeared around the shocks and the rarefaction wave by the first-order time-discretization, the positivity-preserving property is observed in the zoom-in plots. 66 1 0.1 0.8 0.05 0.6 u u 0 0.4 0.2 -0.05 0 -0.1 0 2 4 6 2 2.5 3 3.5 4 x x (a) without limiter (b) zoom in around x = 3 1 0.1 0.8 0.05 0.6 u u 0 0.4 0.2 -0.05 0 -0.1 0 2 4 6 2 2.5 3 3.5 4 x x (c) with limiter (d) zoom in around x = 3 Figure 3.3: Example 3.5.2.3 at T = 0.6. Solid line is the exact solution. Blue squares are point values of the numerical approximation on the Legendre-Gauss- Lobatto points. Mesh with N = 120 cells and CFL = 3.5. 3.5.3 Compressible Euler system At last, we turn to the examination of the applicability of the proposed scheme to the compressible Euler system (3.38). The computational domain Ω = [0, 1] is decomposed into nonuniform mesh. The generic ratio of specific heats is taken to be γ = 1.4. For all the examples, quadratic element is employed with the positivity- CFL preserving limiter on. And the time stepping size is set to be ∆t = kc+|v|k∞ minj hj , 67 where c is the sound speed and v is the velocity. We take CFL = 2 for all the examples. Example 3.5.3.1 (Shock tube). First, let us consider the following shock tube problem.        ρ = 1, if x < 0.5     ρ = 1, if x ≥ 0.5     v = 0, if x < 0.5 ,  v = 0, if x ≥ 0.5         p = 1000,  if x < 0.5 p = 0.01,  if x ≥ 0.5 The numerical approximation is presented in Figure 3.4. The solutions consists of a strong shock wave, a contact discontinuity and a rarefaction wave. All these features are well captured by our scheme. Moreover, with the positivity-preserving limiter on, both the pressure and the density stays positive during the simulation. 1000 6 5 800 4 600 p 3 400 2 200 1 0 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 x x (a) density (b) pressure Figure 3.4: Example 3.5.3.1 at T = 0.01 with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are numerical approximations. Example 3.5.3.2 (Double rarefaction). The double rarefaction problem starts from 68 the following initial condition         ρ = 1, if x < 0.5     ρ = 1, if x ≥ 0.5      v = −2, if x < 0.5 ,  v = 2, if x ≥ 0.5         p = 0.4,  if x < 0.5 p = 0.4,  if x ≥ 0.5 This problem has solution consisting of two symmetric rarefaction waves and trivial contact wave of zero speed. The region between the nonlinear waves around x = 0.5 is close to vacuum, which brings difficulty to the simulation. Without the limiter, the nonlinear solver will experience a hard time to converge. In Figure 3.5, we show the plots for the pressure and the density and we see that the near-vacuum region is well resolved with the help of the positivity-preserving limiter. 0.4 1 0.35 0.8 0.3 0.25 0.6 0.2 p 0.4 0.15 0.1 0.2 0.05 0 0 -0.05 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 x x (a) density (b) Pressure Figure 3.5: Example 3.5.3.2 at T = 0.1, with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are the numerical approximations. Example 3.5.3.3 (Blast wave). The last example is the Sedov point-blast wave [30]. Initially the gas is steady with uniform density one in the whole domain. The pressure is set to be p = 10−9 , except in the central cell, where the pressure is as high as p = 104 . Then a blast-wave starts to propagate from the central cell with 69 a shock front. This problem is difficult, since it involves a low density region and strong shocks. In Figure 3.6 we present the numerical approximation and the exact solution [30] for the pressure and density respectively. The positivity-preserving limiter not only helps keep the low-density region positive, but also add robustness to the nonlinear solver. Without it, the code breaks down due to the failure of convergence of the nonlinear solver. 2000 5 4 1500 3 p 1000 2 500 1 0 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 x x (a) density (b) Pressure Figure 3.6: Example 3.5.3.3 at T = 0.003, with N = 200 and CFL = 2. With the positivity-preserving limiter on. Solid line is the exact solution and blue squares are the numerical approximations. 3.6 Concluding remarks In this chapter, we develop an implicit positivity-preserving DG method with high- order spacial accuracy for one-dimensional conservation laws. This work is an ex- tension of the positivity-preserving limiter in [86, 87] for explicit schemes to implicit ones with backward Euler time discretization. To make the scheme positive, a lower bound for the CFL number is necessary. This conclusion is verified via both theo- retical analysis and numerical experiments. The Crank-Nicolson time discretization 70 is also considered and preliminary analysis on uniform mesh indicates that both the upper bound and lower bound are required to apply the limiter. The positivity- preserving limiter not only makes the numerical approximation physically meaningful but also brings robustness to the scheme and accelerates convergence towards the steady-state solution. The scheme also sees its success on the compressible Euler sys- tem. In the future, we have the following three directions to further explore. First, even though the results in this chapter trivially generalize to multidimensional tensor product meshes and polynomial spaces, positivity-preserving DG schemes for multi- dimensional domain with unstructured mesh needs to be further developed. Second, we plan to generalize the proposed implicit positivity-preserving DG scheme to other types of equations such as the convection-diffusion equations. Thirdly, high-order implicit temporal discretizations need to be considered. Since there are no implicit SSP methods with order greater than one [23], we need to turn to other types of implicit time discretizations, such as BDF methods and fully implicit RK methods. Chapter Four Positivity and Bound-Preserving DG Methods for Relativistic Hydrodynamics 72 4.1 Introduction In this chapter, we discuss discontinuous Galerkin (DG) methods to solve the two- dimensional special relativistic hydrodynamics (RHD), which can be written into a system of conservation laws as below ut + f (u)x + g(u)y = 0, (4.1) with        D   Dw   Dv         m   mw + p   mv  u= , f (u) =  , g(u) =  (4.2)         n   nw   nv + p              E m n as well as its one dimensional version. The method can be easily extended to three- dimensions, but this is not discussed here. In (4.2), p, D, m, n and E are the thermal pressure, mass density, momentum in the x-direction, momentum in the y-direction, and energy, respectively. (w, v) is the velocity field of the fluid. Moreover, units are normalized such that the speed of light is c = 1. If we denote ρ to be the proper rest-mass density, then the conservative variable u can be written as D = γρ, (4.3) m = Dhγw, (4.4) n = Dhγv, (4.5) E = Dhγ − p, (4.6) 73 where γ = (1 − w2 − v 2 )−1/2 is the Lorentz factor and h is the specific enthalpy. To close the system, we specify an equation of state h = h(p, ρ). For ideal gas ρh = ρ + pΓ/(Γ − 1) (4.7) with Γ being the specific heat ratio, such that 1 < Γ ≤ 2 (see for example [66]). Moreover, the sound speed is defined as s r Γp (Γ − 1)(h − 1) cs = = . ρh h Relativistic flows are widely used to model high-energy astrophysical phenomena, such as blast waves of supernova explosions, gravitational collapse and accretion, superluminal jets and gamma-ray bursts. When the speed of the flow is near the speed of the light but there is no strong gravitational field involved, the framework of special relativity is accurate to certain extent to describe the physical phenomena. Physically, the density D and pressure p are positive, and the velocity field (w, v) satisfies w2 + v 2 ≤ 1. Therefore, the admissible set for RHD takes the following form G = {u : D > 0, p(u) > 0, w(u)2 + v(u)2 ≤ 1}. By (4.7), it is easy to see that h > 1, which further yields 0 < cs ≤ 1. It is demonstrated in [42, 77] that G is convex and can be represented in terms of the conservative variables as below √ G = {u : D > 0, E > D2 + m2 + n2 }. (4.8) In order to update the flux in the computation, we further need the inverse of (4.3)- 74 (4.6) from the conservative vector u to the primitive vector w = {ρ, w, v, p}. Unlike its Newtonian counterpart, in RHD we do not have an explicit formula for this inverse map. By (4.6), (4.7) and the definition of γ, we can derive the nonlinear equation satisfied by the pressure p s p m2 + n2 m 2 + n2 f (p; u) := E − −D 1− − = 0, p ∈ [0, +∞) (4.9) Γ−1 (E + p)2 E+p ∂f By a simple calculation one can show that, if u ∈ G, then ∂p (p; u) < 0 and the equation (4.9) has a unique positive solution [77]. In practice, once the bounds (4.8) are satisfied by the numerical solution, (4.9) can be solved efficiently by standard root finding methods. After the pressure is obtained, other quantities can be calcu- lated sequentially and directly via (4.3)-(4.7). Therefore, an efficient and effective conversion from the conservative vector to the primitive vector crucially depends on guaranteeing the bound-preserving property (4.8) for the numerical solution, which is the main objective of this chapter. Numerical simulation of RHD has been intensively studied in the last few decades. The first Eulerian method dates back to the early work by Wilson [74, 75], in which the author used explicit finite differencing techniques and a monotonic transport algorithm to discretize the advection terms of the RHD equations. They applied the von Neumann-Richtmyer artificial viscosity method [71, 60] to handle the shock waves. Though this procedure remained standard through the 1980’s, it turned out to be unable to resolve the extremely strong shock structures that would appear in the ultra-relativistic regime (γ ≥ 2) [7]. A major breakthrough came in 1994 when Mart´ı et al. [35, 36, 37] reformulated the RHD into the conservation form and for the first time introduced the Godunov-type high resolution shock capturing (HRSC) techniques from the classical gas dynamics simulation to that of the RHD. The HRSC techniques produced high-order approximation in the smooth region and were 75 able to capture shocks and steep transients sharply without spurious oscillations. Since then, various Riemann solvers and modern techniques in gas dynamics were extended to the RHD simulation, for example the relativistic Roe solver by Eulderink et al. [20, 21], the HLLE solver extended by Schneider et al. [61] and the recent relativistic HLLC solver carried out by Mignone et al. [42]. Furthermore, Mart´ı and M¨ uller [37] extended the PPM method [14] to the one dimensional RHD and the multidimensional version was accomplished by Mignone et al. [42]. Donat et al. developed a flux splitting method based on the spectral decomposition of the Jacobian matrix in [18]. Kinetic schemes were developed for RHD in [31, 53]. Besides these, other high order methods were also well studied, for example, Dolezal and Wong introduced the essentially non-oscillatory (ENO) scheme for the RHD in 1995 [17], which was extended by Del Zanna and Bucciantini in [83]. Sub- sequently, the WENO algorithm was applied to RHD in [67]. The discontinuous Galerkin methods were also applied to general relativistic hydrodynamics by Radice and Rezzolla [58]. More recently, Zhao and Tang applied WENO limiters to the DG schemes in [93] for the special relativistic case. Moreover, the adaptive mesh refinement (AMR) techniques were proved to be powerful and useful in simulating RHD, for example [84, 73, 25] and many software packages have been developed with AMR support and RHD extension, such as the ENZO [47] and RAMSES [68], etc. Additional methods and details can be found in [80, 81, 76, 92, 34], and the references therein, as well as the review paper [38]. Although the methods mentioned above have been working successfully in most cases, many authors have reported difficulties in maintaining the physical bounds for the numerical approximation, see, e.g. [21, 61, 18, 17, 84, 93], especially in extreme relativistic cases such as flows with large Lorentz factor, strong shocks, low density and pressure. The violation of the physical bounds may lead to the crash 76 √ of the code. Actually, a slight violation of the bound E > D2 + m2 + n2 could √ give negative pressure, and a more severe violation with E < m2 + n2 will even result in the nonexistence of the solution to (4.9), see [61]. Negative pressure ruins the characteristic decomposition in the HRSC methods or kills the Roe solver (see [21]), which all require the square root of the pressure to calculate the speed of the sound. All these could make the code crash in practice. Several ways have been adopted in the literature to get around this problem. For example, in [37], in order to avoid numerical difficulties, the authors set the initial physical quantities, such as the internal energy and pressure, a small number away from zero. In [84], the authors monitor the physical bounds at every time step and once they are broken the calculation will be repeated under a smaller CFL condition with more diffusive schemes. However, these ad hoc techniques are not guaranteed to cure the problem, especially for higher order schemes, and even if they do, the high order accuracy may no longer be maintained. As introduced in Chapter 1, physical bounds preserving high order numerical last few years. For the RHD equation, there are few works to mention. Recently, based on the technique in [26], Wu and Tang [77] developed a bound preserving WENO finite difference scheme for the RHD system (4.1). However, their positivity- preserving technique is based on the flux cut-off limiter proposed by Hu, Adams and Shu in [26]. The main drawback of this method is that to preserve the physical bounds while maintaining the high-order accuracy, a rather severe CFL constraint must be assumed. In this chapter, we restrict ourselves to the DG method and provide a systematic and rigorous way to fix the physical bound violation problem, by extending the bound preserving technique for gas dynamics in [87]. For DG methods, to overcome the oscillations around strong discontinuities, various limiters have been designed, 77 for example the total variation diminishing or total variation bounded (TVD/TVB) limiters in [63, 12] and the WENO limiter in [56]. But for the RHD system in the extreme relativistic cases, the TVD/TVB or the WENO limiter, when applied alone, will not guarantee to preserve the physical bound and the code could still crash. We would like to emphasize that the BP limiter developed in this chapter is not meant as a substitution to the TVB or the WENO limiter or other non- oscillatory techniques, but is simply a remedy to help maintain the physical bounds without affecting the high order accuracy. The TVB or WENO limiter could still be used and may still be necessary in cases involving strong shocks, since the BP limiter is not designed to remove oscillations around shocks, especially when these oscillations are not happening near the physical bounds. However, since the our main objective is the discussion on the bound-preserving limiter, we will not discuss in length about the DG method itself or the TVB or WENO limiters in controlling spurious oscillations. Finally, let us remark that, even though we only discuss the bound-preserving limiter for DG schemes here, the same limiter can be applied to high order finite volume schemes, such as WENO finite volume schemes as well. The organization of this chapter is as follows. In section 4.2, we study the one- dimensional problem, including the details of the limiters, the L1 stability of the method. In sections 4.3 and 4.4, we study the problem in two space dimensions and the implementation of the relativistic axisymmetric jets. Numerical experiments are given in section 4.5. Finally, we will end in section 4.6 with concluding remarks and remarks for future work. 78 4.2 Numerical algorithm in one space dimension In this section, we proceed to construct the bound-preserving DG scheme to solve the one-dimensional relativistic hydrodynamics. 4.2.1 The DG scheme We consider the one-dimensional version of (4.1) on the spatial domain [0, 1] and solve ut + f (u)x = 0, (4.10) where the conservative variable u = (D, m, E)T is defined in (4.3), (4.4), (4.6) with γ = (1 − w2 )−1/2 . The flux is f (u) = (Dw, mw + p, m)T and the equation of state is given in (4.7). The admissible set is defined to be √ G1 = {(D, m, E)T : D > 0, E > D2 + m2 }. It is easy to see that G1 is convex. The DG scheme is defined to seek approximation vector uh (t) ∈ Vh such that in each cell Ij , for any test function v ∈ Vh we have d (uh (t), v)j − (f (uh (t)), vx )j + ˆ − fj+ 1 vj+ ˆ 1 − fj− 1 v + 1 = 0, (4.11) dt 2 2 2 j− 2 − − where vj+ 1 = v(x j+ 1 ), which denotes the left limit of the vector v at xj+ 1 . Likewise 2 2 2 + for vj+ 1. 2 For the numerical flux ˆfj+ 1 = ˆf (uh (x− j+ 1 , t), uh (x+ j+ 1 , t)), we choose the local 2 2 2 79 Lax-Friedrichs flux in the following form ˆ 1 f (a, b) = [f (a) + f (b) − α(b − a)] , (4.12) 2 where α = α(a, b) is a positive real number depending only on the two interface values in the definition of the numerical flux. The value of it is to be chosen by the bound-preserving technique. Other numerical fluxes such as the HLLC Riemann solver [70, 41] could of course also be considered, but will not be discussed here. 4.2.2 Bound-preserving technique For (4.10), direct usage of high order DG methods may result in the appearance of negative density and pressure, and physically irrelevant velocity, leading to ill-posed problems. Moreover, the code may blow up once physically irrelevant quantities appear. Therefore, we would like to apply BP limiters to the scheme. We denote unj and unj to be the numerical solution and its cell average at time level n in cell Ij . For simplicity, throughout the chapter, if we consider generic numerical solution on the whole computational domain Ω, then the subscript j will be omitted. Suppose the exact solution of equation (4.10) is in G1 , we are interested in constructing numerical solutions which are also in G1 . The whole procedure is given below. We start from the forward Euler time discretization, and briefly illustrate the idea. The fully discrete forward Euler DG scheme is defined as follows. In each cell Ij , we seek un+1 h ∈ Vh , such that h i (un+1 , v)j = (unj , v)j ˆn ˆn + λj fj− 1 vj+ 1 − fj+ 1 vj− 1 , ∀v ∈ Vh (4.13) j 2 2 2 2 80 where λj = ∆t/hj and ˆf = ˆf ((unj+ 1 )− , (unj+ 1 )+ ). 2 2 To build the bound-preserving scheme, we follow the general framework outlined in Section 2.2. The first step is to derive sufficient conditions under which given ¯ n+1 unh (x) ∈ G for any x ∈ Ω, we have u j ∈ G for any j. To this end, we generalize the procedure in [87] for the compressible Euler system to the RHD system. First, by taking the test function v = 1 in (4.13) let us write down the equation ¯ n+1 satisfied by the cell average u j , h i ¯ n+1 u = u ¯ n j + λ ˆ f ((u n 1 ) − , (u n 1 ) + ) − ˆ f ((u n 1 )− , (u n 1 )+ ) . (4.14) j j− j− j+ j+ 2 2 2 2 Note that the right hand side involves both the point values u± j± 1 and the cell average 2 ¯ nj . u Then to further get a handle on unh to adjust ¯ n+1 u j , ¯ nj we split the cell average u into summation of point values via Legendre-Gauss-Lobatto rule, which includes the end points xj± 1 . 2 Let ωi be the Legendre Gauss-Lobatto quadrature weights for the interval [− 12 , 12 ] such that M P i=0 ωi = 1, with 2M − 3 ≥ k, and denote the corresponding Gauss- Lobatto points in cell Ij as {xij }M i=1 . Because of the special choice 2M − 3 ≥ k of the Gauss-Lobatto quadrature rule, we have M X ¯ nj = u ωi unj (xji ). i=0 Clearly, unj (xj0 ) = u+ j− 1 and unj (xjM ) = u− j+ 1 . Therefore, by considering ω0 = ωM , we 2 2 have the equation (4.14) can be rewritten into M X   ¯ n+1 u j = ωi unj (xij ) + λj ˆf ((unj− 1 )− , (unj− 1 )+ ) − ˆ f ((unj+ 1 )− , (unj+ 1 )+ ) 2 2 2 2 i=0 81 M −1   X λj ˆ  n − n +  ˆ  n − n +  = ωi unj (xij ) + ω0 (unj− 1 )+ + f (uj− 1 ) , (uj− 1 ) − f (uj− 1 ) , (uj+ 1 ) i=1 ω02 2 2 2 2   n − λj ˆ  n + n −  ˆ  n − n +  + ω0 (uj+ 1 ) + f (uj− 1 ) , (uj+ 1 ) − f (uj+ 1 ) , (uj+ 1 ) 2 ω0 2 2 2 2 M −1   − λj X n i n − n + n = ωi uj (xj ) + ω0 H (uj− 1 ) , (uj− 1 ) , (uj+ 1 ) , i=1 2 2 2 ω0   λj + ω0 H (unj− 1 )+ , (unj+ 1 )− , (unj+ 1 )+ , . (4.15) 2 2 2 ω0 where h i H(a, b, c, η) := b + η ˆf (a, b) − ˆ f (b, c) ¯ n+1 Since unj (xαj ) ∈ G1 and G1 is a convex set, in order to have u j ∈ G1 , we only to require     λj + λj H (unj− 1 )− , (unj− 1 )+ , (unj+ 1 )− , n + n − n , H (uj− 1 ) , (uj+ 1 ) , (uj+ 1 ) , ∈ G1 . 2 2 2 ω0 2 2 2 ω0 (4.16) In general we need to derive condition on η such that given a, b and c ∈ G1 , then we have H(a, b, c, η) ∈ G1 . (4.17) By plugging in the definition of fˆ as in (4.12) into H, we can obtain η H(a, b, c, η) =b + [f (a) + f (b) − α1 (b − a)] 2 η − [f (b) + f (c) − α2 (c − b)] 2 η η η η = [α1 a + f (a)] + (1 − α1 − α2 )b + [α2 c − f (c)] 2 2 2 2 α1 η + η η α2 η − = H (a, α1 ) + (1 − α1 − α2 )b + H (c, α2 ), (4.18) 2 2 2 2 82 where α1 = α(a, b), α2 = α(b, c) and 1 1 H+ (a, α1 ) = a + f (a) , H− (c, α2 ) = c − f (c) α1 α2 Therefore, in order to have (4.17), we need to ensure Condition I Find condition on α such that H ± (u, α) ∈ G1 , for any u ∈ G1 . Condition II Find condition on η such that η η 1 − α1 − α2 ≥ 0 2 2 For the first condition, we have the following lemma. Lemma 4.2.1. Suppose u ∈ G1 and the parameter α satisfies p |w|(h + 1 − 2hτ )γ 2 + τ 4 (h − 1)2 + τ 2 (h − 1)(h + 1 − 2hτ ) α ≥ F (u) = , (4.19) γ 2 (h + 1 − 2hτ ) + τ 2 (h − 1) with τ = (Γ − 1)/Γ, then H± (u, α) ∈ G1 . Proof. For simplicity, we only study the proof for H − , as the proof for H + can be obtained along the same line. Recall that u = (D, m, E)T ∈ G1 and f (u) = (Du, mu + p, m)T . Therefore, αH− (u, α) = αu − f (u) = ((α − u)D, (α − u)m − p, αE − m)T . It is easy to see that to obtain H− ∈ G1 , we need αH− ∈ G1 , i.e. α ≥ c1 := u, (4.20) m α ≥ c2 := , (4.21) E 83 (αE − m)2 ≥ ((α − w)m − p)2 + ((α − w)D)2 . (4.22) Define q(α) = (αE − m)2 − ((α − w)m − p)2 − ((α − w)D)2 . By using (4.3)-(4.6), and define β = α − w, and τ = (Γ − 1)/Γ, we have q(α) = ρ2 (βhγ 2 − ατ (h − 1))2 − β 2 γ 2 − (βhγ 2 w − τ (h − 1))2   = ρ2 (h − 1) β 2 γ 2 (h + 1) + (α2 − 1)τ 2 (h − 1) − 2αβhτ γ 2 + 2βhτ γ 2 w   = ρ2 (h − 1) β 2 γ 2 (h + 1) + (α2 − 1)τ 2 (h − 1) − 2β 2 hτ γ 2   = ρ2 (h − 1) A2 α2 − 2A1 α + A0 ,   where A2 = γ 2 (h + 1 − 2hτ ) + τ 2 (h − 1), (4.23) A1 = uγ 2 (h + 1 − 2hτ ), (4.24) A0 = u2 γ 2 (h + 1 − 2hτ ) − τ 2 (h − 1). (4.25) Since 1 ≤ Γ ≤ 2, then 0 ≤ τ ≤ 1/2. It is easy to check that A2 ≥ 0, q(1) ≥ 0, q(w) ≤ 0, q(−1) ≥ 0. To obtain (4.22), we can take p A1 + (A1 )2 − A2 A0 α ≥ c3 := . (4.26) A2 If u < 0, then c3 ≥ 0 ≥ c1 ≥ c2 . If u ≥ 0, then c3 ≥ c2 ≥ c1 ≥ 0. Therefore, we have H − ∈ G under (4.26). Similarly, we can define A02 = γ 2 (h + 1 − 2hτ ) + τ 2 (h − 1) = A2 , 84 A01 = −uγ 2 (h + 1 − 2hτ ) = −A1 , A00 = u2 γ 2 (h + 1 − 2hτ ) − τ 2 (h − 1) = A0 . to obtain the lower bound of α : p A01 + (A01 )2 − A02 A00 α≥ . (4.27) A02 (4.26) and (4.27) imply p |A1 | + (A1 )2 − A2 A0 α≥ , A2 where A0 , A1 and A2 are given in (4.23)-(4.25). Remark 4.2.1. The same result also holds for the one-dimensional flow with nonzero transverse velocity. The proof is almost the same and is therefore omitted. Remark 4.2.2. In the standard local Lax-Friedrichs flux, the parameter α is cho- sen to be the upper bound for the absolute value of the eigenvalues of the Jacobian ∂f (u)/∂u, which in one dimension is given by (see, e.g. [93]) p p |w| + c 1 − 1/γ 2 + (Γ − 1)(1 − 1/h) α ˜= = p p (4.28) 1 + |w|c 1 + 1 − 1/γ 2 (Γ − 1)(1 − 1/h) In Figure 4.1, we plot α in (4.19) for the bound-preserving requirement (denoted by αbp ) and the spectral radius (4.28) of the Jacobian matrix (denoted by αstandard ), for the two situations with γ = 20 as a constant or with h = 20 as a constant and Γ = 5/3 in both plots. We can see that the wave speed required for the bound- preserving technique is actually slightly lower (that is, less dissipative!) than the standard wave speed used in the Lax-Friedrichs flux. 85 1 0.9995 αbp α αstandard 0.999 0.9985 5 10 15 20 25 30 35 40 h (a) γ=20 1 0.95 αpb 0.9 αstandard α 0.85 0.8 0.75 2 4 6 8 10 γ (b) h=20 Figure 4.1: Plots for the lower bound of α in (4.19) for the bound-preserving re- quirement (denoted by αbp ) and the spectral radius (4.28) of the Jacobian matrix (denoted by αstandard ). The left panel is to keep γ = 20 constant and the right one is to keep h = 20 constant. Given the condition on α in Lemma 4.19, for the condition II, one sufficient condition is to require 1 η≤ (4.29) max{α1 , α2 } 86 In summary, in order to have (4.17) we have the following theorem. Theorem 4.2.1. Suppose a, b and c ∈ G1 , and the parameter α satisfies α1 ≥ max{F (a), F (b)}, α2 ≥ max{F (b), F (c)} with F defined by (4.19), then under the condition (4.29), we have H(a, b, c, η) ∈ G1 Remark 4.2.3. In the very recent work of Wu and Tang [77], which came to our attention after our original submission of [55], there is a similar result as that in Lemma 4.2.1. However, the authors of [77] assume α ≥ α ˜ as defined in (4.28) and show that this condition is sufficient to ensure H± (u, α) ∈ G1 . Our approach appears to be more constructive and provides a less dissipative solver as shown in Figure 4.1. Moreover, for the first order scheme H, the proof in [77] requires a CFL condition 1 η≤ 2 max{α1 ,α2 } to have H ∈ G1 . For our approach, we obtain a more relaxed CFL condition (4.29). Now let us come back to (4.15). In each cell Ij , given unj (x) ∈ G1 for any x ∈ Ij , ¯ n+1 or more relaxed condition unj (xαj ) ∈ G1 for any α, the cell average u j as represented in (4.15) is a convex combination of elements in G1 , since M P α=0 ωα = 1 and ωα > 0. ¯ n+1 And since G1 is a convex set, we must have u j ∈ G1 if (4.16) holds, which by Lemma 4.19 and Theorem 4.2.1 requires 1. The local Lax-Friedrichs parameter satisfies   α (unj+ 1 )− , (unj+ 1 )+ ≥ max{F ((unj+ 1 )− , F ((unj+ 1 )+ )}. 2 2 2 2 87 2. The following CFL type condition ω0 max λj ≤ (4.30) j maxj α((uj+ 1 )− , (unj+ 1 )+ ) n 2 2 Scaling limiter ¯ n+1 We would like to emphasize that the cell average u j is shown to be bound- preserving by the original high order DG scheme, before any limiter is applied. Once the cell average is in control, we can modify the numerical solution through a simple scaling limiter ˜ n+1 ¯ n+1 + θ un+1 ¯ n+1  u j =u j j −u j . (4.31) ˜ n+1 by taking suitable θ ∈ [0, 1], we will have u j ∈ G1 at the Legendre Gauss-Lobatto ˜ n+1 points, and u j is used as the numerical solution at time level n + 1. For scalar equations, we can prove that this modification does not affect the high order accuracy of the original solution un+1 j [88]. Now we give a summary of the complete algorithm. Due to the rounding error, we define     D        √      ε   G1 = u =  : D ≥ ε, E ≥ m 2 + D2 + ε ,  m           E          D        √       ∂Gε1 = u =   m   : D = ε, E = m2 + D2 + ε .         E     then the modification of DG solution unj is given in the following steps. 88 • Set up a small number ε = 10−13 . ¯ jn > ε, then proceed to the following steps. Otherwise, Djn is identified as • If D e nj = unj as the numerical the approximation to vacuum. Therefore, we take u solution and skip the following steps. • Modify density in each Ij : compute bj = minα Djn (xαj ), where {xαj }M α=0 are the Gauss-Lobatto points in the cell Ij . If bj < ε, then take e n = Dn + θD (Dn − Dn ), D j j j j j where n Dj − ε θjD = n . D j − bj e jn as the new numerical density Dn . Then use D j q • Enforce Ejn ˜ n )2 + ε on each Gauss-Lobatto point in each cell ≥ (mnj )2 + (D j Ij : define qαj = unj (xαj ) in cell Ij . If qαj ∈ Gε1 , then take θjα = 1. Otherwise, take θjα to be the root of unj + tqαj = 0,  δ1 (1 − t)¯ t ∈ (0, 1) (4.32) √ where δ1 (u) = E − E 2 + m2 + ε. Then define θj = mini=0,··· ,m θjα , and use e nj = unj + θj (unj − unj ), u as the DG approximation in cell Ij . Remark 4.2.4. In cases where wild data is involved, e.g. the shock heating problem, due to the round-off error, the solution to (4.32) may not strictly make (1 − θij )¯ unj + θjα qαj ∈ Gε1 . In [72], when solving the gas detonation propagation problems with bound 89 preserving DG methods, the authors provide a more robust way to obtain θij without solving any equation. here, we generalize it to the RHD system. It is straightforward to check that for t ∈ (0, 1), unj + tqαj ) ≥ (1 − t)δ1 (¯ δ1 ((1 − t)¯ unj ) + tδ1 (qαj )  given Djn (xαj ) ≥ ε, u ¯ nj ∈ Gε and δ1 qαj < 0, it is sufficient to require unj ) δ1 (¯ (1 − unj ) t)δ1 (¯ + tδ1 (qαj ) = 0, i.e., t= unj ) − δ1 (qαj ) δ1 (¯ unj + tqαj ∈ Gε . Note that the t obtained in this way is in general such that (1 − t)¯ smaller than that by solving the equation (4.32), but the high order accuracy can still be shown following the same way as in [87]. 4.2.2.1 L1 stability Following [79], we can show the L1 -stability of the numerical scheme, for periodic or compactly supported boundary conditions, with the BP limiter. Since Dhn is positive, we have, by the conservative property of the DG scheme, Z Z kDhn kL1 = Dhn (x)dx = Dh0 (x)dx = kDh0 kL1 , Ω Ω where k · kL1 is the standard L1 -norm of w on Ω. Similarly, we can prove kEhn kL1 = kEh0 kL1 . Moreover, it is easy to obtain wΓ m= E, Γ−σ 90 (Γ−1)(h−1) where σ = hγ 2 ≤ Γ − 1. Hence, |w|Γ kmkL1 ≤ kEkL1 ≤ ΓkEkL1 . Γ−σ In the last inequality, we use the fact that |w| ≤ 1. Therefore, kunh kL1 ≤ kDhn kL1 + (Γ + 1)kEhn kL1 = kDh0 kL1 + (Γ + 1)kEh0 kL1 ≤ (Γ + 1)ku0h kL1 , where kukL1 = kDkL1 + kmkL1 + kEkL1 . This implies the L1 -stability of the scheme. 4.2.3 High order time discretizations All the previous analyses are based on first-order Euler forward time discretization. We can also use strong stability preserving (SSP) high-order time discretizations to solve the ODE system ut = Lu. More details of these time discretizations can be found in [65, 64, 23]. In this chapter, we use the third order SSP Runge-Kutta method [65] (1) uh = unh + ∆tL(unh ), (2) 3 1  (1) (1)  uh = unh + uh + ∆tL(uh ) , (4.33) 4 4 1 2  (2) (2)  un+1 = u n + u + ∆tL(u ) , h 3 h 3 h h and the third order SSP multi-step method [65]   n+1 16 n 11 12 u = (u + 3∆tL(un )) + un−3 n−3 + ∆tL(u ) . (4.34) 27 27 11 91 Since an SSP time discretization is a convex combination of Euler forward, by using the limiter designed in section 4.2.2, the numerical solution obtained from the full scheme is also in G1 . 4.3 Numerical algorithm in two space dimensions In this section, we extend the bound-preserving discontinuous Galerkin method to relativistic hydrodynamics in two space dimensions. For simplicity, we use Euler forward time discretization and construct high-order bound-preserving DG schemes to solve (4.1). In this section, we construct the numerical solutions to be in G, which is given in (4.8). For simplicity, we use uniform rectangular meshes. The algorithm including the BP limiter can be easily generalized to unstructured meshes, along the lines h i h i in [91]. The cell is defined as Iij = xi− 1 , xi+ 1 × yj− 1 , yj+ 1 , and the mesh 2 2 2 2 sizes in x and y directions are denoted as ∆x and ∆y, respectively. At time level n, we approximate the exact solution with a vector of polynomials of degree k, n n n unij = (Dij , mnij , nnij , Eijn )T , and define the cell average unij = (Dij , mnij , nnij , E ij )T . Moreover, we denote u+ i− 1 ,j (y), u− i+ 1 ,j (y), u+ i,j− 1 (x), u− i,j+ 1 (x) as the traces of unij on 2 2 2 2 the four edges in cell Iij , respectively. More details can be found in [87]. For simplicity, if we consider a generic numerical solution on the whole computational domain at time level n, then the subscript ij will be omitted. In this section, we only consider high-order schemes, and the one satisfied by the 92 cell averages can be written as Z y 1  ∆t j+ 2    un+1 ij = unij + f u− b 1 (y), u i− 2 ,j + 1 (y) i− 2 ,j − f b u − 1 (y), u i+ 2 ,j + 1 (y) i+ 2 ,j dy ∆x∆y y 1 j− Z x 1 2 ∆t i+ 2    + b u− g i,j− 12 (x), u + i,j− 21 (x) − g b u − i,j+ 12 (x), u + i,j+ 21 (x) dx, (4.35) ∆x∆y x 1 i− 2 f (·, ·) and g where b b(·, ·) are one-dimensional numerical fluxes. For this problem, we still use the one-dimensional local Lax-Friedrichs flux. Suppose (x, y) = (xi− 1 , y0 ) is a point on the vertical cell interface, at which we 2 have two numerical approximations u` = (D` , m` , n` , E` )T and ur = (Dr , mr , nr , Er )T from left and right, respectively. Then the local Lax-Friedrichs flux can be written as 1 f (u` , ur ) = (f (u` ) + f (ur ) − αf (ur − u` )) , b 2 where αf ≥ max{F1 (u` ), F1 (ur )} with F1 is defined by p |w|(h + 1 − 2hτ )γ 2 + τ 4 (h − 1)2 + τ 2 (h − 1)(h + 1 − 2hτ ) F1 (u) = , (4.36) γ 2 (h + 1 − 2hτ ) + τ 2 (h − 1) Γ−1 and the constant τ = Γ . The numerical flux g b can be defined in a similar way with the parameter αg on the horizontal cell interfaces and the corresponding wave speed defined by p |v|(h + 1 − 2hτ )γ 2 + τ 4 (h − 1)2 + τ 2 (h − 1)(h + 1 − 2hτ ) F2 (u) = , (4.37) γ 2 (h + 1 − 2hτ ) + τ 2 (h − 1) We extend the definitions of Hj to two-dimensional problems and define   H1 (u1 , u2 , u3 , λ1 ) = u2 + λ1 bf (u1 , u2 ) − b f (u2 , u3 ) , 93 g (u1 , u2 ) − g H2 (u1 , u2 , u3 , λ2 ) = u2 + λ2 (b b (u2 , u3 )) , ∆t ∆t where λ1 = ∆x and λ2 = ∆y . Following the same proof of Theorem 4.2.1 with some minor changes, we have the following result Lemma 4.3.1. Suppose u1 , u2 , u3 ∈ G, then under a CFL condition max{αf1 , αf2 }λ1 ≤ f (u1 , u2 ) and αf2 is the one for 1, where αf1 is the parameter in the numerical flux b f (u2 , u3 ), then we have b H1 (u1 , u2 , u3 , λ1 ) ∈ G. Similarly, under the CFL condition max{αg1 , αg2 }λ2 ≤ 1 with αg1 and αg2 parameters for g b(u1 , u2 ) and g b(u2 , u3 ), respectively, then we have H2 (u1 , u2 , u3 , λ2 ) ∈ G. To continue, we use L-point Gauss quadratures with L ≥ k + 1 for accuracy to approximate the integrals in (4.35). More details of this requirement can be found h i h i in [12]. The Gauss quadrature points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 are denoted 2 2 2 2 by n o n o β y β pxi = xi : β = 1, · · · , L and pj = yj : β = 1, · · · , L , respectively. Also, we denote wβ as the corresponding weights on the interval  1 1 − 2 , 2 . Different from the notations in previous sections, we use xαi : α = 0, · · · , M } and pˆyj = yˆjα : α = 0, · · · , M pˆxi = {ˆ  h i h i as the Gauss-Lobatto points on xi− 1 , xi+ 1 and yj− 1 , yj+ 1 , respectively. Also, 2 2 2 2 94 we denote να as the corresponding weights on the interval − 21 , 21 .   Then the numerical scheme (4.35) becomes L X h    i un+1 ij = unij + λ1 f u− wβ b 1 i− ,β , u + 1 i− ,β − f b u − 1 i+ ,β , u + 1 i+ ,β 2 2 2 2 β=1 L X h    i − + − + + λ2 b uβ,j− 1 , uβ,j− 1 − g wβ g b uβ,j+ 1 , uβ,j+ 1 , (4.38) 2 2 2 2 β=1 where u− i− 1 ,β = u− (y β ) is a point value in the Gauss quadrature. Likewise for i− 1 ,j j 2 2 the other point values. As the general treatment, we rewrite the cell average on the right hand side as M X X L M X X L unij = να wβ u1αβ = να wβ u2βα , α=0 β=1 α=0 β=1 xαi , yjβ ) and unij (xβi , yˆjα ), respectively. where u1αβ and u2βα denote unij (ˆ i− 21 ,β For each vertical edge, use αf to denote the parameter of the numerical 1 β,j− f (u− flux b i− 1 ,β , u+ i− 1 ,β ) and for the horizontal edge use αg 2 for the parameter of 2 2 b(u− g β,j− 1 , u+ β,j− 1 ). Then define 2 2 i− 21 ,j i− 12 ,β i,j− 12 β,j− 12 αf = max αf , αg = max αg . β=1,··· ,L β=1,··· ,L i− 12 ,j i,j− 21 Let µ = maxi,j αf λ1 + maxi,j αg λ2 , then scheme (4.38) can be written as L M −1   X X − − λ1 un+1 ij = C1 wβ να u1αβ + + ν0 H1 ui− 1 ,β , ui− 1 ,β , ui+ 1 ,β , β=1 α=1 2 2 2 ν0 C 1   λ1 +νM H1 u+ i− 12 ,β , u− i+ 12 ,β , u+ , i+ 12 ,β ν C M 1 L M −1   X X 2 − + − λ2 + C2 wβ να uβα + ν0 H2 uβ,j− 1 , uβ,j− 1 , uβ,j+ 1 , 2 2 2 ν0 C 2 β=1 α=1 95   λ2 +νM H2 u+ β,i− 12 , u− β,i+ 12 , u+ β,j+ 12 , , (4.39) νM C 2 where i− 12 ,j i,j− 12 maxi,j αf λ1 maxi,j αg λ2 C1 = , C2 = . µ µ In (4.39), un+1 ij is the convex combination of u, H1 and H2 . Therefore, we have the following theorem. Theorem 4.3.1. Suppose un ∈ G in scheme (4.38), then un+1 ∈ G, under the CFL condition ∆t i− 1 ,j ∆t i,j− 1 max αf 2 + max αg 2 ≤ ν0 . (4.40) ∆x i,j ∆y i,j Remark 4.3.1. It is straightforward to obtain the bound that F1 (u) ≤ 1, F2 (u) ≤ 1, ∀u ∈ G. i− 12 ,j i,j− 12 In practice, we can replace maxi,j αf and maxi,j αg by 1 to obtain the CFL condition. Based on the above theorem, the numerical cell average we obtain is in G. Of course, the numerical solution un+1 ij might still be placed outside. Hence, we have to modify the numerical solution while keeping the cell average untouched. Due to the rounding error, we define      D                m  √     Gε = u =   : D ≥ ε, E ≥ m2 + n2 + D2 + ε ,      n                   E  96          D             m  √     ε 2 2 2 ∂G = u =   : D = ε, E = D + m + n + ε .      n                   E  Then the modification of unij is given in the following steps. • Set up a small number ε = 10−13 . n n • If Dij > ε, then proceed to the following steps. Otherwise, Dij is identified e nij = unij as the numerical solution as the approximation to vacuum. We take u and skip the following steps. n o n • Modify the density: Compute bij = minαβ Dij xαi , yjβ ), Dij (b n (xβi , ybjα ) . If bij < n ε, then take D e ij as n n D n n D e ij = Dij + θij Dij − Dij , with n D Dij − ε θij = n , Dij − bij n n and use D e ij as the new numerical density Dij . q • Enforce Eijn ≥ ˜ n )2 +  in each cell Iij : Consider u1 and (mnij )2 + (nnij )2 + (D ij αβ u2βα in the cell Iij , respectively. If u1αβ ∈ Gε , then take θαβ 1 = 1. Otherwise, 1 take θαβ to be the root of the equation δ2 (1 − t)unij + tu1αβ = 0,  t ∈ (0, 1). (4.41) √ where δ2 (u) = E − m2 + n2 + D2 + ε. Alternatively, the t can be obtained 2 in the same way as in Remark 4.2.4. Similarly, we can define θβα in the same 97 way for u2βα . Finally, we use e nij = unij + θ(unij − unij ),  1 2 u θ = min θαβ , θβα , α,β as the DG approximation in cell Iij . Now we demonstrate the L1 -stability of the bound-preserving DG scheme. Fol- lowing the same analysis as in Section 4.2, we have kDhn kL1 = kDh0 kL1 and kEhn kL1 = kEh0 kL1 . It is easy to check that wΓ vΓ m= E and n = E, Γ−σ Γ−σ (γ−1)(h−1) where σ = hγ 2 ≤ Γ − 1. Therefore, (|w| + |v|)Γ √ kmkL1 + knkL1 ≤ kEkL1 ≤ 2ΓkEkL1 . Γ−σ In the last inequality, we use the fact that w2 + v 2 ≤ 1. The above analysis yields the L1 -stability of the scheme: √ √ √ kunh kL1 ≤ kDhn kL1 +( 2Γ+1)kEhn kL1 = kDh0 kL1 +( 2Γ+1)kEh0 kL1 ≤ ( 2Γ+1)ku0h kL1 , where kukL1 = kDkL1 + kmkL1 + knkL1 + kEkL1 . 98 4.4 Application to relativistic jets In this section, we study the relativistic axisymmetric jets which is described in the two-dimensional cylindrical coordinates (r, z). We adopt the governing system in [39] and multiply both sides of the system by r. Then it reads ut + f (u)r + g(u)z = s(r, u), (4.42) with          rD   rDu   rDv   0           rm   rmu + rp   rmv   p  u=  , f (u) =   , g(u) =  , s =  .          rn   rnu   rnv + rp   0                  rE rm rn 0 (4.43) We would like to mention that this formulation is not unique. We could also combine the p term in the source and the rp term in the momentum flux to form a non-conservative rpr term, which is closer to the underlying physics. However, since our method is based on the conservative form of the equation we will stay with (4.42). Different from what we have discussed in Section 4.3, the source term s(r, u) is not zero. Therefore, we will demonstrate how to discretize the source terms. The techniques for the discretization of the flux terms and the modification of the numerical approximations follow from the same lines as discussed in Section 4.3, and we will omit them. h i h i In this section, the cell is defined as Iij = ri− 1 , ri+ 1 × zj− 1 , zj+ 1 , with 2 2 2 2 99 1 ≤ i ≤ Nr , 1 ≤ j ≤ Nz and the mesh sizes in r and z directions are denoted as ∆r and ∆z, respectively. If not otherwise stated in this section, the notations follow those in Section 4.3. We consider high order schemes only, the one satisfied by the cell averages can be written as Z z 1  1 n ∆t j+ 2    un+1 ij = uij + f u− b 1 (z), u i− 2 ,j + 1 (z) i− 2 ,j − f b u − 1 (z), u i+ 2 ,j + 1 (z) i+ 2 ,j dz 2 ∆r∆z z 1 j− Z r 1  2 ∆t i+ 2    + b u− g i,j− 12 (r), u + i,j− 21 (r) − g b u − i,j+ 12 (r), u + i,j+ 12 (r) dr, ∆r∆z r 1 i− 2 1 n + uij + ∆tsnij (4.44) 2 where snij is the cell average of s in the cell Iij at time level n. We use L-point Gauss quadratures with L ≥ k+1 2 to approximate snij . The Gauss quadrature points h i h i on ri− 1 , ri+ 1 and zj− 1 , zj+ 1 are denoted by 2 2 2 2 n o n o pri = riβ : β = 1, · · · , L and pzj = zjβ : β = 1, · · · , L , respectively. Also, we denote wβ as the corresponding weights on the interval  1 1 − 2 , 2 . Then L X L 1 n X n u + ∆tsij = wα wβ Hs (uαβ αβ ij , sij , ∆t), 2 ij α=1 β=1 where uαβ n α β αβ n α αβ ij = uij (ri , zj ), sij = sij (ri , uij ) and Hs (u, s, ∆t) = 1 2 u + ∆ts. Since s is a function of u and r, Hs can be written as a function of u, r and ∆t, i.e. H e s (u, r, ∆t) = Hs (u, s, ∆t). We can choose ∆t sufficiently small, such that e s (u, r, ∆t) = Hs (u, s, ∆t) ∈ G, and the result is given in the following lemma. H Lemma 4.4.1. Suppose u ∈ G, if the time mesh size ∆t satisfies r ∆t ≤ αe(u, r), (4.45) 2 100 where p (wΓhγ 2 )2 − (Γ − 1)2 (h − 1)2 − Γ2 γ 2 (h − 1)2 + 2Γhγ 2 (h − 1) − uhΓγ 2 α e(u, r) = , (Γ − 1)(h − 1) e s (u, r, ∆t) ∈ G. then H Proof. It is easy to see that e s (w, r, ∆t) = (rD, rm + ∆tp, rn, rE)T . H Define q(w, r, ∆t) = (rE)2 − (rD)2 − (rm + ∆tp)2 − (rn)2 . Since r > 0, we only need to find ∆t such that q(w, r, ∆t) ≥ 0. It is not difficult to check that   Γ E= −1 p (4.46) σ uΓ m= p (4.47) σ vΓ n= p (4.48) σ γΓ D= p, (4.49) (Γ − 1)(h − 1) (Γ−1)(h−1) where σ = hγ 2 . Therefore, we have q(w, r, ∆t) = (rE)2 − (rD)2 − (rm + ∆tp)2 − (rn)2 " 2  2  2  2 # 2 rΓ rγΓ uΓr rΓv =p −r − − + ∆t − σ (Γ − 1)(h − 1) σ σ " # 2 Γ2 2Γ ∆t2 2uΓ ∆t   2 2 Γγ =r p +1− − 2 − − σ2γ 2 σ r σ r (Γ − 1)(h − 1)  2  ∆t ∆t = −r2 p2 2 + 2A +B , r r 101 where uΓhγ 2 A= (Γ − 1)(h − 1) and Γhγ 2 (Γ − 2) − Γ2 γ 2 B= − 1 < −1 (Γ − 1)2 (h − 1) Let q(w, r, ∆t) ≥ 0, we need √ ∆t ≤ ( A2 − B − A)r. (4.50) With the above lemma, we can state the following theorem Theorem 4.4.1. Suppose un ∈ G, then un+1 ∈ G under the conditions riα ∆t ≤ max e(uαβ α α ij , ri ), 2 1 ≤ α, β ≤ L 1 ≤ i ≤ Nr 1 ≤ j ≤ Nz and ∆t i− 1 ,j ∆t i,j− 1 ν0 max αf 2 + max αg 2 ≤ . ∆r i,j ∆z i,j 2 i− 12 ,j i,j− 12 where αf and αg have the same definition as those in (4.40). Now we have finished all the theoretical analysis. 102 4.5 Numerical experiments In this section, we present numerical examples in both one and two dimensions to verify the bound-preserving property of the proposed method. In order to demon- strate the effectiveness of the proposed BP limiter, we have deliberately not applied any other non-oscillatory limiters, such as the TVD/TVB or WENO limiters. As expected from the theory, the numerical solution stays in the physical bounds and computations could proceed stably, even though there are spurious oscillations in some test results due to the lack of non-oscillatory limiters. We should emphasize that this bound-preserving limiter, which is extremely local (implemented completely inside each cell without using information from neighboring cells) and inexpensive, is not meant to substitute other non-oscillatory limiters or techniques. It can be used together with other non-oscillatory limiters or techniques, such as TVD/TVB or WENO limiters, or even WENO finite volume schemes, to obtain both bound- preserving and non-oscillatory performance, but we will not pursue such approach here in order to be focused on the bound-preserving property. We would also like to emphasize that many of the examples are chosen to contain solutions particularly challenging in terms of the bound-preserving requirement, such that even with the usual TVD/TVB or WENO limiters added, the code would still fail without the proposed bound-preserving technique, due to the breaking of the physical bound which causes ill-posedness of the problem as well as the nonexistence of the solution to (4.9). Unless otherwise indicated, we consider the third-order RKDG method (k = 2) with the local Lax-Friedrichs numerical flux (4.12). The CFL number is set to be 0.15. The specific heat ratio Γ is taken to be 5/3, unless otherwise stated. 103 4.5.1 One-dimensional experiments Example 4.5.1.1 (Smooth flow). Consider a one-dimensional flow in the domain Ω = [0, 1] with the following initial states ρ0 (x) = 1 + 0.9999999 sin(2πx), u0 (x) = 0.9, p0 (x) = 1.0 The boundary condition is set to be periodic. The exact solution is ρ(x, t) = 1 + 0.9999999 sin(2π(x − 0.9t)), w(x, t) = 0.9, p(x, t) = 1.0 In Table 4.1, we present the numerical results for the proposed method with and without the bound-preserving limiter. For the case k = 3, we take ∆t = CFL(∆x)4/3 so that the time error will not dominate. Since this is a high speed smooth flow with its lowest density near zero, the bound-preserving limiter does get turned on. As observed and discussed in [90] and [79], the SSP Runge-Kutta (RK) methods might degenerate the accuracy when the bound-preserving limiter is applied, while for the SSP multi-step (M-S) method, the full order of accuracy is recovered. Moreover, in Figure 4.2, we plot the CPU time against the L2 error for polynomial degrees k = 1, 2, 3 with and without the BP limiter. The following two observations could be made: (i) As one of the advantages of high order DG methods, for a given error, the computational time spent is diminishing with increasing order; (ii) The additional cost introduced by the BP limiter is very small. Example 4.5.1.2 (Moderate blast wave). The relativistic blast wave problems are standard tests (see [38]) for a numerical relativistic hydrodynamical code. The com- putational domain is Ω = [0, 1]. We first consider a moderate case, which has the 104 no limiter with limiter (RK) with limiter (M-S) k h L2 error order L2 error order L2 error order 1/20 3.50 e-2 – 4.85 e-2 – 4.55 e-2 – 1/40 8.72 e-3 2.00 1.06 e-2 2.11 1.05 e-2 2.12 1 1/80 2.17 e-3 2.01 2.55 e-3 2.05 2.55 e-3 2.03 1/160 5.43 e-4 2.00 6.13 e-4 2.03 6.11 e-4 2.06 1/320 1.36 e-4 2.00 1.50 e-4 2.04 1.49 e-4 2.03 1/20 2.16 e-3 – 3.47 e-3 – 2.40 e-3 – 1/40 2.77 e-4 2.97 5.05 e-4 2.78 2.79 e-4 3.10 2 1/80 3.48 e-5 2.99 9.42 e-5 2.42 3.48 e-5 3.00 1/160 4.35 e-6 3.00 1.87 e-5 2.33 4.35 e-6 3.00 1/320 5.44 e-7 3.00 3.68 e-6 2.35 5.44 e-7 3.00 1/20 8.99 e-5 – 1.82 e-4 – 1.52 e-4 – 1/40 5.60 e-6 4.00 2.09 e-5 3.12 5.82 e-6 8.03 3 1/80 3.46 e-7 4.02 2.30 e-6 3.18 3.50 e-7 4.05 1/160 2.16 e-8 4.00 2.70 e-7 3.09 2.19 e-8 4.00 1/320 1.35 e-9 4.00 1.86 e-8 3.86 1.46 e-9 3.91 Table 4.1: Example 4.5.1.1: One-dimensional accuracy test at T = 0.4 for the second- , third- and fourth-order DG methods with and without the BP limiters. The CFL number for the multistep method is one-third of that for the RK method. k is the degree of polynomial and h is the meshsize. initial states   (10.0, 0.0, 13.33), x < 0.5  (ρ0 , u0 , p0 ) = (4.51)  (1.0, 0.0, 0.0),  x > 0.5 The initial discontinuity gives rise to a transonic rarefaction wave propagating left, a shock wave propagating right and a contact discontinuity in between. The numerical approximation in Figure 4.3 shows that our method can resolve these structures quite well. Small oscillations are observed around the contact discontinuity, since we have not applied any non-oscillatory limiters such as the TVD/TVB or WENO limiters. For the same reason, in some of the following examples, we also observe spurious oscillations due to the lack of non-oscillatory limiters. Example 4.5.1.3 (Strong relativistic blast wave). Next we examine a much more 105 k=1 k=1, with limiter k=2 10-2 k=2, with limiter k=3 -3 k=3, with limiter 10 -4 10 L2 error -5 10 10-6 10-7 -8 10 -9 10 10-1 100 101 102 CPU time Figure 4.2: Example 4.5.1.1: Log-log plot of the CPU time against the L2 error for polynomial degrees k = 1, 2, 3 with and without the BP limiter. challenging blast wave example, of which the initial states are   (1.0, 0.0, 103 ), x < 0.5  (ρ0 , u0 , p0 ) = (4.52)  (1.0, 0.0, 0.0), x > 0.5  This example is more relativistic than the previous one. The high relativity is due to the large enthalpy of the left state, which is h ' 2.5 × 103  1. This results in a thermodynamically relativistic configuration. The structure of the solution is the same as the moderate case, except for the formation of a very thin dense shell behind the shock in the density and a highly curved profile for the rarefaction fan in the velocity. The relativistic shock propagating at a Lorentz factor γ ' 6 [84]. In Figure 4.4, despite of small oscillations, we can see that our method resolves the curved part in the profile of the velocity very well and captures the thin shell in the density with little smearing. To further test the bound preserving property of our method, we consider a more 106 10 14 12 8 10 6 8 ρ p 6 4 4 2 2 0 0 0.2 0.4 0.6 0.8 1 x 0 0.2 0.4 0.6 0.8 1 x (a) density ρ (b) pressure p 0.8 0.7 0.6 0.5 0.4 u 0.3 0.2 0.1 0 -0.1 0 0.2 0.4 0.6 0.8 1 x (c) velocity w Figure 4.3: Example 4.5.1.2: Blast wave with initial condition (4.51) at T = 0.5 on the mesh of 200 cells, approximated by the third-order RKDG method with the BP limiter. The squares represent the approximate cell averages and the solid line is the exact solution. extreme example. The initial condition now is   (1.0, 0.0, 104 ), x < 0.5  (ρ0 , u0 , p0 ) = (4.53)  (1.0, 0.0, 0.0), x > 0.5  with the specific heat ratio Γ = 4/3. In this case the enthalpy of the left state is h ' 4 × 104 . In the profile of the density, the width of the thin shell is approximately 107 1000 10 8 800 6 600 p ρ 4 400 2 200 0 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 x x (a) density ρ (b) pressure p 1 0.8 0.6 u 0.4 0.2 0 0.2 0.4 0.6 0.8 x (c) velocity w Figure 4.4: Example 4.5.1.3: Strong blast wave with initial condition (4.52) at T = 0.4 approximated by the third-order RKDG method with the BP limiter on the mesh with 200 cells. The square represents the approximate cell averages and the solid line is the exact solution. 2.5 × 10−3 and the rarefaction part in the velocity becomes even more curved. In Figure 4.5, we can see the good performance of our proposed method. Example 4.5.1.4 (Shock-heating problem). The shock-heating problem (see [83, 84, 37, 38]) is a standard benchmark problem to examine the ability of the method to deal with strong shocks. A warm wall is located at x = 0 and cold gas flows in at x = 1.0. When the gas gets reflected on the wall, a reverse strong shock forms 108 15 15 10 10 ρ ρ 5 5 0 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 x x (a) density ρ, 400 cells (b) density ρ, 800 cells 1 1 0.8 0.8 0.6 0.6 u u 0.4 0.4 0.2 0.2 0 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 x x (c) velocity w, 400 cells (d) velocity w , 800 cells Figure 4.5: Example 4.5.1.3: Strong blast wave with initial condition (4.53) at T = 0.3 approximated by the third-order RKDG method with the BP limiter. The mesh is decomposed into 400 cells (left column) and 800 cells (right column). The square represents the approximate cell averages and the solid line is the exact solution. and propagates to the right. In our experiment, the inflow velocity is set to be vin = 0.999999, with the corresponding Lorentz factor γ = 707.1. The initial density is ρ0 = 1.0 and the initial pressure p0 = 0.0. The specific heat ratio in this example is Γ = 4/3. The analytical solution can be found in [37]. In Figure 4.6, we observe that, due to the BP limiter, near the physical bounds, there is no oscillation and the bounds are well preserved. The shocks are sharply 109 captured and the oscillations behind the shock may be eliminated with the non- oscillatory techniques such as the TVB or the WENO limiter. There is an undershoot in the left end point of the density, which is the so-called “wall-heating” effect. All these show the effectiveness of the proposed BP limiter in the ultra-relativistic regime with very strong shocks. 3 7 6 2.5 5 2 p/100000 4 ρ/1000 1.5 3 1 2 0.5 1 0 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8 x x (a) density ρ (b) pressure p 1 0.8 0.6 u 0.4 0.2 0 0.2 0.4 0.6 0.8 x (c) velocity w Figure 4.6: Example 4.5.1.4: One-dimensional shock-heating problem on the mesh of 200 cells, at T = 1.5 with vin = 0.999999, approximated by the third-order RKDG method with the BP limiter. The squares represent the numerical results and the solid line is the exact solution. Example 4.5.1.5 (One-dimensional Riemann problem with non-zero transverse ve- locity). The last one-dimensional test is the Riemann problem with non-zero trans- 110 verse velocity, see [84] and [42]. The domain is still Ω = [0, 1] and the initial states are   (1.0, 0.0, 0.9, 103 ), x < 0.5  (ρ0 , u0 , v0 , p0 ) = (4.54)  (1.0, 0.0, 0.9, 10−2 ), x > 0.5  Compared with Example 4.5.1.3, where there is no transverse velocity, the density has a smaller jump and a wider dense shell, however the distance between the tail of the rarefaction and the contact discontinuity is much smaller, which requires very high resolution to resolve the structure. The purpose of this example is to show the robustness of the BP limiter even when the transverse speed is nonzero and close to the speed of light. Admittedly, to correctly restore the right wave speed, it is unavoidable to refine the mesh (see the comparison between 400 cells and 6400 cells in Figures 4.7 and 4.8 and also results in [84] with the adaptive mesh refinement technique), but this goes beyond the out purpose and we will not pursue it here. Our results on meshes of 400 and 6400 cells are presented in Figure 4.7 and Figure 4.8. These results match those in [84] and [42] well. We further test this problem with a transverse velocity v = 0.999. In this case, the dense shell in the density is much thinner and the transverse velocity increases from v = 0.999 to v = 0.99967 at the rarefaction part, which corresponds to a Lorentz factor γ ' 39. We see in Figure 4.9, our method is still robust in this severe case. 111 5 5 4 4 3 3 ρ ρ 2 2 1 1 0 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 x x (a) density ρ, 400 cells (b) density ρ, 6400 cells 1 0.95 0.95 0.9 0.9 v 0.85 0.85 v 0.8 0.8 0.75 0.75 0.7 0.7 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 x x (c) velocity v, 400 cells (d) velocity v, 6400 cells Figure 4.7: Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity at T = 0.6, approximated by the third-order RKDG method with the BP limiter. Approximations of the density and the transverse velocity on meshes of 400 and 6400 cells. The squares represent the approximate cell averages and the solid lines are the exact solutions. 4.5.2 Two-dimensional experiments Example 4.5.2.1 (Smooth flow). First, the accuracy of the method in 2-D is checked √ by considering the following smooth flow in the domain Ω = [0, 2]2 . ρ = 1 + 0.999999 sin[2π(r − vt) · b], vx = 0.9, vy = 0.2, p = 1.0 (4.55) 112 1000 1000 800 800 600 600 p p 400 400 200 200 0 0 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 x x (a) pressure p, 400 cells (b) pressure p, 6400 cells 0.4 0.35 0.3 0.3 0.25 0.2 0.2 u u 0.15 0.1 0.1 0.05 0 0 -0.05 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 x x (c) velocity w, 400 cells (d) velocity w, 6400 cells Figure 4.8: Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity at T = 0.6, approximated by the third-order RKDG method with the BP limiter. The approximations of w and p on meshes of 400 and 6400 cells are presented respectively. The squares represent the approximate cell averages and the solid lines are the exact solutions. where r = (x, y), v = (vx , vy ) and b = (cos α, sin α) is the direction vector, along which the wave propagates. In this example we take α = π/4. In Table 4.2, the numerical errors in L2 and L∞ norms of the RKDG method with the BP limiter are listed. We can see that, in this case, the limiter preserves the high-order accuracy. Example 4.5.2.2 (Oblique shock wave). Next, let us check a two-dimensional √ oblique 1D shock tube test in the domain [0, 2/2]2 . We consider the strong blast 113 5 1000 4 800 3 600 p ρ 2 400 1 200 0 0 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 x x (a) density ρ (b) pressure p 1 0.03 0.9995 0.999 0.02 0.9985 u v 0.998 0.01 0.9975 0.997 0 0.9965 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 x x (c) velocity w (d) velocity v Figure 4.9: Example 4.5.1.5: One-dimensional Riemann problem with non-zero transverse velocity v = 0.999 at T = 0.6, approximated by the third-order RKDG method with the BP limiter. The approximations of the w, v, p and ρ on 3200 cells are presented. The blue solid lines are the numerical results and the black ones are the exact solutions. wave in Example 4.5.1.3 propagating along the direction at a 45 degree angle. Ini- tially, the state is divided into two parts by the line x + y = 1 and is set to be   (1.0, 0.0, 0.0, 103 ), x + y < 1  (ρ0 , u0 , v0 , p0 ) = . (4.56)  (1.0, 0.0, 0.0, 0.0), x + y > 1  The domain is divided into 120 × 120 cells and the terminal time T = 0.3. In Figure 114 k h L2 error order L∞ error order 1/10 4.53 e-2 – 1.98 e-1 – 1/20 2.61 e-3 4.11 1.64 e-2 3.60 2 1/40 2.97 e-4 3.14 1.65 e-3 3.31 1/80 3.84 e-5 2.95 2.16 e-4 2.93 1/160 5.13 e-6 2.91 2.75 e-5 2.97 Table 4.2: Example 4.5.1.1: Two-dimensional accuracy test at T = 0.2 for the third- order RKDG method with the BP limiter. k is the degree of polynomial and h is the mesh size. 4.10, we present the contours as well as the cuts along the x = y line of the numerical √ approximation for ρ, p and w2 + v 2 . We see that as for the 1D counterpart, despite of small oscillations, all the structures: the shock, contact continuity and rarefaction, are well resolved and the physical bounds are preserved. Example 4.5.2.3 (Two-dimensional Riemann problem). The two dimensional Rie- mann problem involves the interactions of elementary waves, which initially separate four constant states. In [83], the authors extended the 2-D Riemann problem from the classical Newtonian hydrodynamics [62] to the relativistic flows. We take the same initial condition as in [83].      (0.1, 0.0, 0.0, 10−2 ), x, y > 0    (0.1, 0.99, 0.0, 1.0), x < 0 < y  (ρ0 , u0 , v0 , p0 ) =     (0.5, 0.0, 0.0, 1.0), x, y < 0     (0.1, 0.0, 0.99, 1.0), y<0