Monotonicity-preserving finite element methods for hyperbolic problems
Abstract
This thesis covers the development of monotonicity preserving finite element methods for hyperbolic problems. In particular, scalar convection-diffusion and Euler equations are used as model problems for the discussion in this dissertation. A novel artificial diffusion stabilization method has been proposed for scalar problems. This technique is proved to yield monotonic solutions, to be \ac{led}, Lipschitz continuous, and linearity preserving. These properties are satisfied in multiple dimensions and for general meshes. However, these results are limited to first order Lagrangian finite elements. A modification of this stabilization operator that is twice differentiable has been also proposed. With this regularized operator, nonlinear convergence is notably improved, while the stability properties remain unaltered (at least, in a weak sense). An extension of this stabilization method to high-order discretizations has also been proposed. In particular, arbitrary order space-time isogeometric analysis is used for this purpose. It has been proved that this scheme yields solutions that satisfy a global space-time discrete maximum principle unconditionally. A partitioned approach has also been proposed. This strategy reduces the computational cost of the scheme, while it preserves all stability properties. A regularization of this stabilization operator has also been developed. As for the first order finite element method, it improves the nonlinear convergence without harming the stability properties. An extension to Euler equations has also been pursued. In this case, instead of monotonicity-preserving, the developed scheme is local bounds preserving. Following the previous works, a regularized differentiable version has also been proposed. In addition, a continuation method using the parameters introduced for the regularization has been used. In this case, not only the nonlinear convergence is improved, but also the robustness of the method. However, the improvement in nonlinear convergence is limited to moderate tolerances and it is not as notable as for the scalar problem. Finally, the stabilized schemes proposed had been adapted to adaptive mesh refinement discretizations. In particular, nonconforming hierarchical octree-based meshes have been used. Using these settings, the efficiency of solving a monotonicity-preserving high-order stiff nonlinear problem has been assessed. Given a specific accuracy, the computational time required for solving the high-order problem is compared to the one required for solving a low-order problem (easy to converge) in a much finer adapted mesh. In addition, an error estimator based on the stabilization terms has been proposed and tested. The performance of all proposed schemes has been assessed using several numerical tests and solving various benchmark problems. The obtained results have been commented and included in the dissertation.
Full text
Monotonicity-preserving finite element methods for hyperbolic problems Jesús Bonilla ADVERTIMENT La consulta d’aquesta tesi queda condicionada a l’acceptació de les següents condicions d'ús: La difusió d’aquesta tesi per mitjà del r e p o s i t o r i i n s t i t u c i o n a l UPCommons (http://upcommons.upc.edu/tesis) i el repositori cooperatiu TDX ( h t t p : / / w w w . t d x . c a t / ) ha estat autoritzada pels titulars dels drets de propietat intel·lectual únicament per a usos privats emmarcats en activitats d’investigació i docència. No s’autoritza la seva reproducció amb finalitats de lucre ni la seva difusió i posada a disposició des d’un lloc aliè al servei UPCommons o TDX. No s’autoritza la presentació del seu contingut en una finestra o marc aliè a UPCommons (framing). Aquesta reserva de drets afecta tant al resum de presentació de la tesi com als seus continguts. En la utilització o cita de parts de la tesi és obligat indicar el nom de la persona autora. ADVERTENCIA La consulta de esta tesis queda condicionada a la aceptación de las siguientes condiciones de uso: La difusión de esta tesis por medio del repositorio institucional UPCommons (http://upcommons.upc.edu/tesis) y el repositorio cooperativo TDR (http://www.tdx.cat/?localeattribute=es) ha sido autorizada por los titulares de los derechos de propiedad intelectual únicamente para usos privados enmarcados en actividades de investigación y docencia. No se autoriza su reproducción con finalidades de lucro ni su difusión y puesta a disposición desde un sitio ajeno al servicio UPCommons No se autoriza la presentación de su contenido en una ventana o marco ajeno a UPCommons (framing). Esta reserva de derechos afecta tanto al resumen de presentación de la tesis como a sus contenidos. En la utilización o cita de partes de la tesis es obligado indicar el nombre de la persona autora. WARNING On having consulted this thesis you’re accepting the following use conditions: Spreading this thesis by the i n s t i t u t i o n a l r e p o s i t o r y UPCommons (http://upcommons.upc.edu/tesis) and the cooperative repository TDX (http://www.tdx.cat/?localeattribute=en) has been authorized by the titular of the intellectual property rights only for private uses placed in investigation and teaching activities. Reproduction with lucrative aims is not authorized neither its spreading nor availability from a site foreign to the UPCommons service. Introducing its content in a window or frame foreign to the UPCommons service is not authorized (framing). These rights affect to the presentation summary of the thesis as well as to its contents. In the using or citation of parts of the thesis it’s obliged to indicate the name of the author.
Universitat Politècnica de Catalunya Doctoral Thesis Monotonicity-preserving finite element methods for hyperbolic problems Author: Jesús Bonilla Supervisor: Prof. Santiago Badia A thesis submitted in fulfillment of the requirements for the degree of Doctor of Philosophy in the Doctorat en Enginyeria Civil Departament d’Enginyeria Civil i Ambiental Barcelona, November 2019
To my family & friends iii
“A thesis has to be presentable... but don’t attach too much importance to it. If you do succeed in the sciences, you will do later on better things and then it will be of little moment. If you don’t succeed in the sciences, it doesn’t matter at all. ” Paul Ehrenfest v
Abstract This thesis covers the development of monotonicity-preserving finite element methods for hyperbolic problems. In particular, scalar convection-diffusion and Euler equations are used as model problems for the discussion in this dissertation. A novel artificial diffusion stabilization method has been proposed for scalar problems. This technique is proved to yield monotonic solutions, to be local extremum diminishing (LED), Lipschitz continuous, and linearity preserving. These properties are satisfied in multiple dimensions and for general meshes. However, these results are limited to first order Lagrangian finite elements. A modification of this stabilization operator that is twice differentiable has been also proposed. With this regularized operator, nonlinear convergence is notably improved, while the stability properties remain unaltered (at least, in a weak sense). An extension of this stabilization method to high-order discretizations has also been proposed. In particular, arbitrary order space-time isogeometric analysis is used for this purpose. It has been proved that this scheme yields solutions that satisfy a global space-time discrete maximum principle unconditionally. A partitioned approach has also been proposed. This strategy reduces the computational cost of the scheme, while it preserves all stability properties. A regularization of this stabilization operator has also been developed. As for the first order finite element method, it improves the nonlinear convergence without harming the stability properties. An extension to Euler equations has also been pursued. In this case, instead of monotonicity-preserving, the developed scheme is local bounds preserving. Following the previous works, a regularized differentiable version has also been proposed. In addition, a continuation method using the parameters introduced for the regularization has been used. In this case, not only the nonlinear convergence is improved, but also the robustness of the method. However, the improvement in nonlinear convergence is limited to moderate tolerances and it is not as notable as for the scalar problem. Finally, the stabilized schemes proposed had been adapted to adaptive mesh refinement discretizations. In particular, nonconforming hierarchical octree-based meshes have been used. Using these settings, the efficiency of solving a monotonicity-preserving high-order stiff nonlinear problem has been assessed. Given a specific accuracy, the computational time required for solving the high-order problem is compared to the one required for solving a low-order problem (easy to converge) in a much finer adapted mesh. In addition, an error estimator based on the stabilization terms has been proposed and tested. The performance of all proposed schemes has been assessed using several numerical tests and solving various benchmark problems. The obtained results have been commented and included in the dissertation. vii
5.3 Nonlinear stabilization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 5.3.1 Differentiable stabilization . . . . . . . . . . . . . . . . . . . . . . . 103 5.4 Adaptive mesh refinement . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 5.4.1 Error estimators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.4.2 Refinement strategy . . . . . . . . . . . . . . . . . . . . . . . . . . 107 5.5 Nonlinearsolver.................................107 5.6 Numericalresults................................109 5.6.1 Convergence...............................109 5.6.2 Linear discontinuity . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.6.3 Circular discontinuity . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.6.4 Compression corner . . . . . . . . . . . . . . . . . . . . . . . . . . 115 5.6.5 Reflectedshock.............................117 5.7 Conclusions...................................120 6 Conclusions and future work 121 6.1 Conclusions...................................121 6.2 Futurework...................................123 Bibliography 125 xv
List of Figures 2.1 Representation of the symmetric node jsym of jwith respect to i. . . . . . 15 2.2 Convergence test, L2(Ω) error versus size of the mesh. For P1and Q1FE meshes ranging from h= 1/12 to h= 1/96. Newton’s method has been used with parameters q= 4,ε= 10−7,σ=|v|h410−8and γ= 10−10. . . . 29 2.3 Stabilized solution of the straight propagation of a discontinuity test using the steady version of discrete problem (2.12) with two stabilization choices (2.26)or(2.7). ................................. 30 2.4 Stabilized solution of the straight propagation of a discontinuity test using the steady version of discrete problem (2.12) with two stabilization choices (2.26) and (2.7). The stabilization parameters used for the smoothed version are q= 25,ε= 10−4,σ=|v|10−9, and γ= 10−10. ......... 30 2.5 Straight propagation test solution at the outflow boundary ∂Ω\Γin. Using the steady version of discrete problem (2.12) and nonlinear diffusion (2.24), for different values of qand ε,σ=|v|ε10−5,γ= 10−10, and both nonlinear solvers in Sect. 2.8. The result in brackets shows the number of iterations if no projection to Vadm hisdone................... 31 2.6 Evolution of global discrete maximum principle (DMP) violation during nonlinear iterations when avoiding the projection step in Algs. 1 and 2 for the straight propagation of a discontinuity test. . . . . . . . . . . . . . 32 2.7 Straight propagation test nonlinear iterations as mesh refined from 12 × 12 Q1to 96 ×96 Q1, for both Alg. 1 and Alg. 2. The shock capturing parameters used are q= 4,ε= 10−2,σ=|v|h410−6, and γ= 10−10. . . . 33 2.8 Circular propagation test solution at the outflow boundary ∂Ω\Γin. Using the steady version of discrete problem (2.12) and nonlinear diffusion (2.24), for different values of qand ε,σ=|v|ε10−5,γ= 10−10 and both nonlinear solvers in Sect. 2.8. The result in brackets shows the number of iterations if no projection to Vadm hisdone................... 34 2.9 Stabilized solution of the circular convection test using the steady version of the discrete problem (2.12) and the nonlinear diffusion (2.24) for two different parameter choices. . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.10 Evolution of global DMP violation during nonlinear iterations when avoiding the projection step in Algs. 1 and 2 for the circular propagation of a discontinuity. Using q= 25,ε= 10−4,σ=|v|10−9,γ= 10−10. . . . . . . . 35 xvii
2.11 Cross-sections of each for the figures rotated in the three body rotation benchmark. The parameters used are q= 25,γ= 10−8,σ=|v|10−10, ε= 10−4, and ∆t= 10−3, in a 150 ×150 Q1element mesh. The discrete problem (2.12) is used in combination with three different artificial diffusions (2.24) and (2.6) leading to a LED scheme, and (2.15) leading to a globalDMPscheme. .............................. 36 2.12 3 Body rotation test results using discrete problem (2.12) and two different artificial diffusions ((2.24) leading an LED scheme, and (2.15) with (2.26) leading a global DMP scheme). Using a 150 ×150 Q1element mesh, and parameters: q= 25,γ= 10−8,σ=|v|10−10,ε= 10−4, and ∆t= 10−3. . . 37 2.13 Burger’s equation solutions at t= 0.5using discrete problem (2.12) and (2.6) with (2.24). Using a 150 ×150 Q1element mesh, ∆t= 10−2, and two sets of parameters q,γ,σ, and ε...................... 38 3.1 Representation of the basis functions of V2 hin one dimension, with its associated Greville abscissae. . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Representation of the polytope Qiin two dimensions, the symmetric node xsym ij of xjwith respect to xi,xaand xb. .................. 47 3.3 Second and third order discretizations obtained from the k-refinement of an initial first order discretization. Notice that shape functions are depicted for interior knots, at boundary knots shape functions become interpolatory, see Fig. 3.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.4 Convergence in time results for problem (3.11), using standard and partitioned space-time schemes. . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.5 Convergence in space results for problem (3.12). . . . . . . . . . . . . . . . 55 3.6 Effect of the regularization parameters for first order discretizations. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picardscheme. ................................. 56 3.7 Effect of the regularization parameters for second order discretization. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picardscheme. ................................. 57 3.8 Solution of problem (3.13) at t= 0.5for first to fourth order discretizations. 58 3.9 Solution of problem (3.13) at t= 0.5for first to fourth order discretizations. 58 3.10 Solution of problem (3.14) using scheme (3.6), and different discretization orders....................................... 59 xviii
3.11 Three body rotation test initial conditions. . . . . . . . . . . . . . . . . . 60 3.12 Three body rotation test results at t= 1 using scheme (3.6), q= 10, σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 200×200×1000 control points is used with 250 subdomains in the temporal direction. .................................... 61 3.13 Three body rotation test results at t= 1 using scheme (3.6), q= 10, σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 125 subdomains in the temporal direction havebeenused.................................. 61 3.14 Three body rotation test results at t= 1 using scheme (3.6), q= 10, σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 250 subdomains in the temporal direction havebeenused.................................. 62 3.15 Three body rotation test profiles for t= 1 at p(x−0.5)2+ (y−0.5)2= 0.25 using scheme (3.6), q= 10,σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 125 and 250 subdomains in the temporal direction have been used. . . . . . . . . . . . 62 4.1 usym drawing................................... 73 4.2 Compression corner scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.3 Density convergence for successive mesh refinements. . . . . . . . . . . . . 80 4.4 Reflected shock scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.5 Reflected shock convergence history for q= 1................. 82 4.6 Reflected shock convergence history for q= 2................. 82 4.8 Reflected shock convergence history for q= 10. ............... 82 4.7 Reflected shock convergence history for q= 5................. 83 4.9 Sod shock initial condition and solution for the differentiable scheme using parameters q= 10,σ= 10−3,ε= 10−5, and γ= 10−10............ 83 4.10 Comparison of L1error and computational cost (total number of iterations) for different regularization parameters choices at the Sod’s shock test. ....................................... 84 4.11 Scramjet test scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.12 Scramjet Mach contours when a mesh of 63695 Q1elements is used, with parameters q= 5,γ= 10−10, and ˜ε= 1. ................... 86 4.13 Scramjet Mach contours when a mesh of 63695 Q1elements is used, with parameters q= 5,γ= 10−10, and ˜ε= 1. ................... 86 4.15 Scramjet Mach contours when a mesh of 18476 Q1elements is used, with parameters q= 2,γ= 10−10, and ˜ε= 1. ................... 86 xix
4.14 Scramjet Mach contours when a mesh of 63695 Q1elements is used, with parameters q= 2,γ= 10−10, and ˜ε= 1. ................... 87 4.16 Comparison of the convergence behavior for the Scramjet test and different regularization parameters choices. A coarse mesh of 18476 Q1elements is used........................................ 87 4.17 Comparison of the convergence behavior for the Scramjet test and different regularization parameters choices. A fine mesh of 63695 Q1elements is used........................................ 88 5.1 Example of a mesh with hanging nodes. ................... 94 5.2 usym drawing ..................................100 5.3 Compression corner scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . 109 5.4 Convergence of ku−uhkL1(Ω) to a solution with a discontinuity. . . . . . . 110 5.5 Evolution of the mesh refinement process. ˜ηKwith high-order scheme is used in the left column. Low-order scheme with Kelly estimator is used in the central column. ˜ηKwith low-order scheme is used in the right column. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 9, and 9. For the high-order with q= 10, the refinement steps are 1, 2, 3, 5, and 5. . . . . . . . . . . . . . . . . . . . . . 111 5.6 Time and elements convergence comparison for the transport problem with a linear discontinuity, q= 1........................112 5.7 Time and elements convergence comparison for the transport problem with a linear discontinuity, q= 2........................112 5.8 Time and elements convergence comparison for the transport problem with a linear discontinuity, q= 10. ......................113 5.9 Evolution of the mesh refinement process. ˜ηKwith high-order scheme is used in the left column. Low-order scheme with Kelly estimator is used in the central column. ˜ηKwith low-order scheme is used in the right column. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 7, and 7. For the high-order with q= 10, the refinement steps are 1, 2, 3, 4, and 4. . . . . . . . . . . . . . . . . . . . . . 114 5.10 Time and elements convergence comparison for the transport problem with a circular convection field, q= 1. ....................115 5.11 Time and elements convergence comparison for the transport problem with a circular convection field, q= 2. ....................115 5.12 Time and elements convergence comparison for the transport problem with a circular convection field, q= 10.....................115 5.13 Evolution of the mesh refinement process. ˜ηKwith high-order (right) and low-order (left) schemes are used. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 8, and 8. For the high-order with q= 10, the refinement steps are 1, 2, 3, 4, and 4. . 116 xx
5.14 Time and elements convergence comparison for the compression corner problem......................................117 5.15 Reflected shock scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.16 Time and elements convergence comparison for the reflected shock problem.118 5.17 Evolution of the mesh refinement process. ˜ηKwith low-order scheme is used. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 4, 5, 6, and 7. The lower two figures are the high-order with q= 2 (top) and low-order (bottom) results at their last refinement step. . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 xxi
List of Tables 2.1 Straight propagation test errors and iterations, using the steady version of discrete problem (2.12) and nonlinear diffusion (2.24), for different values of qand ε,σ=|v|ε10−5,γ= 10−10, and both nonlinear solvers in Sect. 2.8......................................... 31 2.2 Circular propagation test errors and iterations, using the steady version of discrete problem (2.12) and nonlinear diffusion (2.24), for different values of qand ε,σ=|v|ε10−5,γ= 10−10, and both nonlinear solvers in Sect. 2.8......................................... 34 3.1 Measured convergence rates in L2norm and H1semi-norm, for problem (3.11)....................................... 54 3.2 Measured convergence rates in L2and H1norms, for problem (3.12). . . . 55 4.1 Experimental convergence rates for both problems. . . . . . . . . . . . . . 81 4.2 Reflected shock solution values at every region. . . . . . . . . . . . . . . . 81 4.3 Domain coordinates for the scramjet test. . . . . . . . . . . . . . . . . . . 85 5.1 Reflected shock solution values at every region. . . . . . . . . . . . . . . . 117 xxiii
1.2. Thesis objectives 5 cannot benefit from the higher accuracy of high-order FEs, which is of special interest for problems that combine shocks and smooth regions. Moreover, in the context of hp-adaptive schemes, forcing p= 1 in the vicinity of discontinuities and shocks might become cumbersome. Moreover, many authors refer to strong stability preserving (SSP) Runge Kutta (RK) methods to achieve high-order convergence in time. However, to achieve high convergence rates these kind of methods require to satisfy a CFL-like condition [37]. The motivation of using an implicit time integrator was precisely to avoid stability conditions on the time step length. Therefore, we will explore alternatives to SSP RK methods to achieve high-order time integration. •Extension to first order hyperbolic systems of equations. After exploring the previously mentioned goals, an important step is to start working with systems of equations. In particular, to extend at least some of the previous achievements to a problem closer to the application in the motivation. As a first step, we will consider the extension of the methods developed for scalar problems to Euler equations. •Extension to AMR FE schemes. In the case of problems with discontinuities, the ability of automatically adapt the resolution of the mesh to the features of the problem can notably improve the convergence. Therefore, we consider important to ensure that the methods developed in this thesis are well suited to be used in this kind of discretizations. Moreover, the stabilization methods explored in this thesis are characterized by restricting its action to the vicinity of discontinuities. Hence, we will explore the possibility of using this property in the AMR process to identify which regions of the mesh need to be refined. •Assessment of the efficiency of high-order monotonicity-preserving schemes in AMR context. As previously mentioned, using this kind of methods requires to solve a stiff nonlinear problem. Previous goals explicitly attempt to improve the nonlinear convergence. However, we find interesting to check whether it is still clearly better to use a high-order scheme with AMR. That is, for a given accuracy, we would like to test whether it is more efficient to use a high-order method (with a stiff nonlinear problem), or if it is better to use a low-order method with a much finner mesh. •Code development. All the results in this thesis (but the ones in Chapter 2) have been implemented in the in-house FE library FEMPAR [9, 11]. FEMPAR stands for Finite Element Multiphysics PARallel solvers, and it is a finite element library that provides all the
6Chapter 1. Introduction necessary tools for computing FE approximations of PDEs, e.g. from discretization to numerical linear algebra. It is important to mention that FEMPAR is a collaborative software project and without the contributions of current and former developers the results in this thesis would have been impossible to achieve. 1.3 Document structure The first chapter of this thesis contains a brief introduction, the motivations of the research developed, and the specific goals of this thesis. Chapters 2 to 5 contain the main contributions of this study. Each one of this chapters corresponds directly to the publications in the list of the next section. The chapters are self-contained, preserve the structure of the paper, and can be read independently. However, we have tried to keep the notation as homogeneous as possible. Chapter 2 is devoted to the development of a monotonicity-preserving stabilization method for first order Lagrangian FEs in arbitrary meshes. It also contains a differentiable version that improves the computational cost. First, Sect. 2.1 contains an introduction to monotonicity-preserving stabilization methods for FEs. In Sect. 2.2, the continuous problem and its discretization using the FE method are presented. Sect. 2.3 contains the formulation of the novel nonlinear stabilization method. Sect. 2.4 is devoted to the monotonicity analysis of the proposed method. An alternative approach is presented in Sect. 2.5. Lipschitz continuity of the methods is proved in Sect. 2.6. A differentiable version the previous method is presented in Sect. 2.7. Sect. 2.8 is devoted to nonlinear solvers. Different numerical experiments are introduced in Sect. 2.9. Finally, we draw some conclusions in Sect. 2.10. Chapter 3 contains an extension of the previous stabilization method to space–time arbitrary high-order isogeometric analysis. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems in Sect. 3.2. Then, the stabilization techniques are introduced in Sect. 3.3. Sect. 3.4 is devoted to a partitioned time integration scheme. Afterwards, we introduce a regularized version of the stabilization term in Sect. 3.5. Finally, we show numerical experiments in Sect. 3.6 and draw some concluding remarks in Sect. 3.7. In Chapter 4 the previous differentiable shock detector techniques are combined with stabilization methods for systems of equations. The resulting scheme is proved to be local bounds preserving for first order hyperbolic problems. In Sect. 4.2 we present the continuous Galerkin (cG) discretization for scalar convection and Euler equations. Sect. 4.3 is devoted to the definition of the stabilization terms. We describe the nonlinear solvers used in Sect. 4.4. Then, we present the numerical experiments performed in Sect. 4.5. Finally, we draw some conclusions in Sect. 4.6. In Chapter 5 the performance of the methods developed in Chapter 2 and Chapter 4 are evaluated in the context of AMR. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems and hyperbolic systems in Sect. 5.2.
1.4. List of publications and conference participations 7 Then, the stabilization techniques are introduced in Sect. 5.3. Sect. 5.4 is devoted to the AMR strategy. Afterwards, we introduce the nonlinear solvers used in Sect. 5.5. Finally, we show numerical experiments in Sect. 5.6 and draw some conclusions in Sect. 5.7. Finally, Chapter 6 summarizes the conclusions and the main goals achieved in the thesis at hand. In addition, we also introduce possible future works to pursue based on the developments in this thesis. 1.4 List of publications and conference participations The work presented in this thesis has also been published in international peer reviewed journals, as well as in international conferences. The journal articles written in the scope of this thesis are listed below. [4] S. Badia and J. Bonilla,Monotonicity-preserving finite element schemes based on differentiable nonlinear stabilization, Computer Methods in Applied Mechanics and Engineering 313 (2017) 133–158. [16] J. Bonilla and S.Badia,Maximum principle preserving space-time isogeometric analysis, Computer Methods in Applied Mechanics and Engineering 354 (2019) 422–440. [6] S. Badia, J. Bonilla, S. Mabuza and J. Shadid,Differentiable local bounds preserving stabilization for first order hyperbolic problems. Submitted. [17] J. Bonilla and S. Badia,Monotonicity-preserving finite element schemes with adaptive mesh refinement for hyperbolic problems. In preparation. In addition, the following article was also written during this thesis. However, for the forthcoming chapters we will focus only on cG methods. [5] S. Badia J. Bonilla, and A. Hierro,Differentiable monotonicity-preserving schemes for discontinuous Galerkin methods on arbitrary meshes, Computer Methods in Applied Mechanics and Engineering 320 (2017) 582–605. Moreover the developments described in this thesis have been presented in the following international conferences. 2016 S. Badia∗and J. Bonilla,Monotonicity preserving nonlinear stabilization for hyperbolic scalar problems, Conference on the Mathematics of Finite Elements and Applications. Uxbridge, England. 2017 S. Badia∗and J. Bonilla,Finite element methods preserving maximum principles, Finite Elements in Fluids Conference. Rome, Italy. 2017 S. Badia and J. Bonilla∗,Monotonicity preserving finite element methods for scalar convection–diffusion problems, European Workshop on High Order Nonlinear Numerical Methods for Evolutionary PDEs. Stuttgart, Germany.
8Chapter 1. Introduction 2017 J. Bonilla∗and S. Badia,High-order monotonicity preserving finite element methods for scalar convection–diffusion problems, European Conference on Numerical Mathematics and Advanced Applications. Voss, Norway. 2019 J. Bonilla∗and S. Badia,Monotonicity preserving stabilization for convection dominated flows, International Congress on Industrial and Applied Mathematics. Valencia, Spain It is worth mentioning that during the course of this doctoral studies the author performed a six-months research stay at Sandia National Laboratories. Among other tasks, the author could collaborate with Prof. Shadid and his team, which led to the following proceeding as well as the journal article [6]. 2018 J. Bonilla, S. Mabuza, J.N. Shadid, and S. Badia,On differentiable linearity and local bounds preserving stabilization methods for first order conservation law systems, (2018).
Chapter 2 Monotonicity preserving stabilization for linear FEs This chapter is focused on a nonlinear stabilization technique for scalar conservation laws with implicit time stepping. The method relies on an artificial diffusion method, based on a graph-Laplacian operator. It is nonlinear, since it depends on a shock detector. Further, the resulting method is linearity preserving. The same shock detector is used to gradually lump the mass matrix. The resulting method is LED, positivity preserving, and also satisfies a global DMP. Lipschitz continuity has also been proved. However, the resulting scheme is highly nonlinear, leading to very poor nonlinear convergence rates. We propose a smooth version of the scheme, which leads to twice differentiable nonlinear stabilization schemes. It allows one to straightforwardly use Newton’s method and obtain quadratic convergence. In the numerical experiments, steady and transient linear transport, and transient Burgers’ equation have been considered in 2D. Using the Newton method with a smooth version of the scheme we can reduce 10 to 20 times the number of iterations of Anderson acceleration with the original non-smooth scheme. In any case, these properties are only true for the converged solution, but not for iterates. In this sense, we have also proposed the concept of projected nonlinear solvers, where a projection step is performed at the end of every nonlinear iterations onto a FE space of admissible solutions. The space of admissible solutions is the one that satisfies the desired monotonic properties (maximum principle or positivity). 2.1 Introduction Many PDEs satisfy some sort of maximum principle or positivity property. However, numerical discretizations usually violate these structural properties at the discrete level, with implications in terms of accuracy and stability, e.g., leading to non-physical local oscillations. It is well-understood now how to build methods that satisfy some sort of DMP based on explicit time integration combined with FVM or discontinuous Galerkin (dG) schemes [27, 68]. However, implicit time integration is preferred in problems with multiple scales in time when the fastest scales are not relevant. E.g., under-resolved simulations of 9
10 Chapter 2. Monotonicity preserving stabilization for linear FEs multi-scale problems in time are essential in plasma physics [54]. Unfortunately, implicit DMP-preserving hyperbolic solvers are scarce and not so well developed. In the frame of FE discretizations, the local instabilities present in the solution of hyperbolic problems have motivated the use of so-called shock capturing schemes based on artificial diffusion (see, e.g., [49]). These methods introduce nonlinear stabilization, in contrast with classical SUPG-type linear stabilization techniques [47, 48]. Since linear schemes are at most first-order accurate and highly dissipative [36], recent research on FE techniques for conservation laws has focused on the development of less dissipative nonlinear schemes. Many of these ideas come from the numerical approximation of convection dominated CDR, where one encounters similar issues. The cornerstone of these methods is the design of a nonlinear artificial diffusion that vanishes in smooth regions and works on discontinuities or sharp layers. Many residual-based diffusion methods have been considered so far (see, e.g., [32] and references therein). Most of these approaches have failed to reach DMP-preserving methods. A salient exception is the method by Burman and Ern [25], which satisfies a DMP under mesh restrictions. Recently, due to some interesting novel approaches in the field, the state-of-the-art in nonlinear stabilization has certainly advanced [7, 13, 14, 24, 26, 64, 65]. Implicit FE schemes for hyperbolic problems rely on four key ingredients: 1. The first ingredient is the definition of the shock detector that only activates the nonlinear diffusion around shocks/discontinuities. Recent nonlinear stabilization techniques have been developed based on shock detectors driven by gradient jumps [7, 23] or edge differences [13, 64, 65]. The use of such schemes was proposed in [23] for 1D problems and extended to multiple dimensions in [7]. A salient property of the scheme in [7] is that it is DMP-preserving, but it relies on the DMP of the Poisson operator, which is only true under stringent constraints on the mesh. Another salient feature of the gradient-jump diffusion approach in [7] is the fact that it leads to so-called linearity preserving methods, i.e., the artificial diffusion vanishes for first order polynomials. This property is related to high-order convergence on smooth regions [66]. A modification of the nonlinear diffusion in [64] that also satisfies this property is proposed in [65]. 2. The second ingredient is the amount of diffusion to be introduced on shocks, which is the amount of diffusion introduced in a first order linear scheme. In this sense, one can consider flux-corrected transport techniques [67]. 3. The third ingredient is the form of the discrete viscous operator. In order to keep the DMP on arbitrary meshes, Guermond and Nazarov have proposed to use graph-theoretic, instead of PDE-based, operators for the artificial diffusion terms. This approach has been used in [13, 93] (for the steady-state convection–diffusion– reaction problem) and in [39] (for linear conservation laws) combined with artificial diffusion definitions similar to the one in [38].
2.1. Introduction 11 4. The fourth ingredient is the perturbation of the mass matrix, in order to satisfy a DMP. Full mass lumping is one choice, but it introduces an unacceptable phase error. For continuous FE methods, improved techniques can be found in [40]. Alternatively, limiting-type strategies are used, e.g., in [64, 65]. 5. The method in [13] is Lipschitz continuous, which is needed for the well-posedness of the resulting nonlinear scheme. However, in practice, all the methods presented above are still highly nonlinear, and nonlinear convergence becomes very hard and expensive. It leads to a fifth additional ingredient that has not been considered so far in much detail. In order to reduce the computational cost of these schemes, we consider the smoothing of the nonlinear artificial diffusion, to make it differentiable up to some fixed order. The possibility to define smooth nonlinear schemes can improve the nonlinear convergence of the methods and make them practical for realistic applications. Further, the smoothing step enables advanced linearization strategies based on Newton’s method. It also involves the development of efficient nonlinear solvers, e.g., based on the combination of Newton, line search, and/or Anderson acceleration techniques. All the results commented above are restricted to linear (or bilinear) FEs. We are not aware of the existence of high-order implicit DMP-preserving FE schemes. For explicit time integration and limiters, second order methods can be found in [39]. The use of hp-adaptive schemes that keep first order schemes around shocks has been proposed in [44]. In this chapter, we propose a novel nonlinear stabilization method that satisfies a DMP, positivity, and LED properties at the discrete level. It combines: (1) a novel shock detector related to the one in [7], which is simple and linearity preserving; (2) the graph- Laplacian artificial viscous term proposed in [38]; (3) an edge FCT-type definition of the amount of diffusion (see [64]); (4) a novel gradual mass lumping technique that exploits the same shock detector used for the artificial diffusion. We prove that the resulting method ticks all the boxes, i.e., it is total variation diminishing (TVD), DMP, positivitypreserving, linearity preserving, Lipschitz continuous, and introduces low dissipation. With regard to the last point, we prove that the amount of diffusion is the minimum needed in our analysis to prove the DMP. Further, we consider a novel approach to design a smoothed version of the resulting scheme that is twice differentiable. We prove that linear preservation is weakly enforced in this case, but all the other properties remain unchanged. Finally, we analyze the effect of the smoothing in the computational cost, and observe a clear reduction in the CPU cost of the nonlinear solver when using the smooth version of the method proposed herein while keeping almost unchanged the sharp layers of the non-smooth version. Future work will be focused on the entropy stability analysis of these schemes for nonlinear scalar conservation laws. A partial result in this direction is the proof of entropy stability for a related method when applied to the 1D Burger’s equations (see [23]).
12 Chapter 2. Monotonicity preserving stabilization for linear FEs This chapter is structured as follows. In Sect. 2.2 the continuous problem and its discretization using the FE method are presented. Sect. 2.3 contains the formulation of a novel nonlinear stabilization method. Sect. 2.4 is devoted to the monotonicity analysis of the proposed method. An alternative approach is presented in Sect. 2.5. Lipschitz continuity of the methods is proved in Sect. 2.6. A differentiable version the previous method is presented in Sect. 2.7. Sect. 2.8 is devoted to nonlinear solvers. Different numerical experiments are introduced in Sect. 2.9. Finally, in Sect. 2.10 we draw some conclusions. 2.2 Preliminaries 2.2.1 The continuous problem Let Ω⊂Rdbe a bounded domain, where dis the space dimension, and (0, T]the time interval. The scalar conservation equation reads: find u(x, t)such that ∂tu+∇·f(u) = g, on Ω×(0, T], where f∈Lip(R;Rd)is the flux. It is also subject to the initial condition u(x,0) = u0∈L∞(Ω) and boundary condition u(x, t) = u(x, t)on the inflow Γin . ={(x, t)∈∂Ω× (0, T]|f(x, t)·n(x, t)<0}. There exist a unique entropy solution uof the above problem that satisfies the entropy inequalities ∂tE(u) + ∇·F(u)≤0for all convex entropies E∈Lip(R;R)with its associated entropy fluxes Fi(u) = Ru 0E0(v)f0 i(v)dv, 1≤i≤d (see Kružkov [55]). Let us consider the weak form of this problem consists in seeking u such that u=uon Γin ×(0, T]and (∂tu, v)+(f0(u)·∇u, v) = (g, v)∀v∈L2(Ω),(2.1) almost everywhere in (0, T], with g∈L2(Ω). 2.2.2 Finite element spaces and meshes Let Thbe a conforming partition of Ωinto elements, K. Elements can be triangles or quadrilaterals for d= 2, or tetrahedrals or hexahedra for d= 3. The set of interpolation nodes of This represented by Nh, whereas Nh(K)denotes the set of nodes belonging to element K∈ Th. Moreover, Ωiis the macroelement composed by the union of the elements K∈ Thsuch that i∈ Nh(K).Nh(Ωi)denotes the set of nodes in that macroelement. The continuous linear FE space is defined as Vh. =vh∈ C0(Ω) : vh|K∈Pk(K)∀K∈ Th for triangular or tetrahedral elements (replacing P1(K)by Q1(K)for quadrilateral or hexahedral elements). P1(K)(resp., Q1(K)) is the space of polynomials with total (resp.,
2.2. Preliminaries 13 partial) degree less or equal to 1. The nodal basis of Vhis written {ϕi}i∈Nh, and the FE functions can be expressed as vh=Pi∈Nhϕivi, where viis the value of vhat node i. 2.2.3 The semi-discrete problem The semi-discrete Galerkin FE approximation of (2.1) reads: find uh∈Vhsuch that uh(Γin, t) = πh(u)and (∂tuh, vh)+(f0(uh)·∇uh, vh)=(g, vh)∀vh∈Vh,(2.2) for t∈(0, T], with initial conditions uh(·,0) = πh(u0).πhdenotes a FE interpolation, e.g., the Scott-Zhang projector [83]. Using the notation Muh. = (uh,·)and F(wh)uh. = (f0(wh)·∇uh,·)we can write problem (2.2) in compact form as M∂tuh+F(uh)uh=g(2.3) in V0 h, i.e., the dual space of Vh. Further, we define Mij . = (ϕj, ϕi),Fij(uh). = (f0(uh)· ∇ϕj, ϕi), and gi. = (g, ϕi). In order to carry out the time discretization of (2.3), let us consider a partition of the time domain (0, T]into sub-intervals (tn, tn+1], with 0. =t0< t1< . . . < tN. =T. We consider the Backward-Euler (BE) implicit time integrator to keep at the time-discrete level the monotonicity properties of the semi-discrete problem, leading to the discrete problem: given u0 h . =πh(u0)∈Vh, compute for n= 1, . . . , N −1 Mδtun+1 h+F(un+1 h)un+1 h=gin V0 h,(2.4) where δtun+1 h . = ∆t−1 n+1(un+1 h−un h), and ∆tn+1 . =|tn+1 −tn|. Implicit strong stability preserving Runge-Kutta methods [53] also preserve the monotonic properties at the discrete level [53], under some restrictions on the time step size. For the sake of brevity we consider the BE scheme. Systems (2.3) and (2.4) will be supplemented with additional stabilization terms to minimize the oscillations generated by the Galerkin FE approximation. Of particular interest are methods which provide solutions that satisfy the following property for all nodes, for zero forcing terms. Definition 2.2.1 (Local DMP).A solution u∈Vhsatisfies the local DMP if umin i≤ui≤umax i,where umax i . = max j∈Nh(Ωi)\{i}uj, umin i . = min j∈Nh(Ωi)\{i}uj. Actually, for steady problems, if this is satisfied for all i∈ Nh, then the extrema will be at the boundary and there exist no local extrema. Furthermore, it is useful to define local extremum diminishing (LED) methods for transient problems.
14 Chapter 2. Monotonicity preserving stabilization for linear FEs Definition 2.2.2 (LED).A method is called LED if for g= 0 and any time in (0, T], the solution satisfies dtui≤0if uiis a maximum and dtui≥0if uiis a minimum. For time-discrete methods, the same definition applies, replacing dtby δt. 2.3 Nonlinear stabilization We want to design a linearity preserving LED method for stabilizing the scalar semidiscrete hyperbolic problem (2.3) (or the discrete problem (2.4)), described in the previous section. As written above, this method is based on a graph-theoretic approach. Let us consider a nonlinear stabilization operator B(uh) : Vh→V0 hand denote Bij(uh). = hB(uh)ϕj, ϕii. Particularly, we require that the stabilization term will satisfy the following properties (see also [38]): 1. compact support: Bij(uh)=0if j /∈ Nh(Ωi)for any uh∈Vh, 2. symmetry: Bij(uh) = Bji(uh)for any uh∈Vh, 3. conservation: Pj6=iBij(uh) = −Bii(uh)for any uh∈Vh, 4. linear preservation: B(uh) = 0 for any uh∈P1(Ω). To achieve this properties we define the nonlinear stabilization term hB(wh)uh, vhi. =X i∈NhX j∈Nh(Ωi) νij(wh)viuj`(i, j), uh, vh∈Vh,(2.5) where the graph-theoretic Laplacian is defined as `(i, j). = 2δij −1, and the artificial diffusion computed as νij(wh). = max {αi(wh)Fij(wh), αj(wh)Fji(wh),0}for i6=j, νii(wh). =X j∈Nh(Ωi) j6=i νij(wh),(2.6) where αi(·)is the shock detector. We note that this choice leads to a symmetric stabilization operator B(wh). In order to define the shock detector, let us introduce some notation. Let i∈ Nhbe a node of the mesh, va vector field, and wa scalar field. Let rij =xj−xibe the vector pointing from nodes ito jin Nhand ˆ rij . =rij |rij |. Let xsym ij be the point at the intersection between the line that passes through xiand xjand ∂Ωi that is not xj(see Fig. 2.1). The set of all symmetric nodes with respect to node iis represented with Nsym h(Ωi). We define rsym ij . =xsym ij −xi, and usym j . =uh(xsym ij ). Then, one can define the jump and the mean of the unknown gradient at node iin direction
2.6. Lipschitz continuity 21 2.6 Lipschitz continuity In the next, we want to prove the Lipschitz continuity of the nonlinear operator at every time step, i.e., T:Vh→V0 hdefined as T(uh). = ∆t−1 n+1M(uh)uh+K(uh)uh−g−∆t−1 n+1M(uh)un h. In order to prove the Lipschitz continuity of T(·), we must deal with the nonlinear stabilization and gradual mass lumping terms. The Galerkin terms can be handled using the fact that f∈Lip(R;Rd). Let us introduce the following semi-norm generated by the graph-Laplacian operator |w|`. =v u u t 1 2X i∈NhX j∈Nh(Ωi) (wi−wj)2. Further, we define |v|as the supremum of |f(v)|for v∈Vadm h, where Vadm h⊂Vhis the subspace of functions that satisfy the global DMP in Def. 2.5.1. Theorem 2.6.1. Let us consider a non-degenerate partition Th. Given un h∈Vhand g∈V0 h, the nonlinear operators B(·) : Vh→V0 hand M(·) : Vh→V0 hare Lipschitz continuous in Vadm hfor q∈N+, since they satisfy hB(u)−B(v), zi ≤ qhd−1|v||u−v|`|z|`,for any z∈Vh, hM(u)−M(v), zi ≤ C(qhd 2|u−v|`+ku−vk)kzk,for any z∈Vh. Proof. We assume that the FE mesh is quasi-uniform in order to reduce technicalities. However, the proof for Lipschitz continuity can be extended to more general meshes. We denote A=cB as AhBand A < cB as A.B, for any positive constant cthat does not depend on the numerical or physical parameters. From the definition of the nonlinear stabilization in (2.5), we get |hB(u)u, zi−hB(v)v, zi| ≤X i∈NhX j∈Nh(Ωi) νij(v)`(i, j)(uj−vj)zi (2.16) +X i∈NhX j∈Nh(Ωi) (νij(u)−νij(v))`(i, j)ujzi . Using the definition of |v|, the Cauchy-Schwarz inequality, the fact that kϕik ≤ Chd/2, and the inverse inequality k∇vhk.h−1kvhkfor vh∈Vh(see [21]), we get: Fij(w)≤ |v|k∇ϕik2 Lkϕjk2 L.hd−1|v|,(2.17)
22 Chapter 2. Monotonicity preserving stabilization for linear FEs for any w∈Vadm h. Using (2.17), the first term in the RHS of (2.16) is bounded as follows: X i∈NhX j∈Nh(Ωi) νij(v)`(i, j)(uj−vj)zi .hd−1|v||u−v|`|z|`. The second term is bounded using the Cauchy-Schwarz inequality: X i∈NhX j∈Nh(Ωi) (νij(u)−νij(v))`(i, j)ujzi(2.18) .X i∈NhX j∈Nh(Ωi) 1 2(νij(u)−νij(v))2(ui−uj)2 1 2 ×|z|`. Using (2.17), we have: νij(u)−νij(v) = max{αi(u)Fij(u), αj(u)Fji(u),0}−max{αi(v)Fij(v), αj(v)Fji(v),0}(2.19) ≤max{(αi(u)Fij(u)−αi(v)Fij(v), αj(u)Fji(u)−αj(v)Fji(v),0} .hd−1|v|max{|αi(u)−αi(v)|,|αj(u)−αj(v)|}. Let us assume that Pj∈Nh(Ωi){{|∇uh·rij|}}ij 6= 0. (The other case is straightforward.) On one hand, for a non-degenerate FE mesh, we have that ch ≤rij ≤Ch,j∈ Nsym h(Ωi), for positive constants c, C that do not depend on h. Using this fact in the definition of the shock detector (2.7), we get: αi(u)1 q=Pj∈Nh(Ωi)J∇uhKij Pj∈Nh(Ωi)2{{|∇uh·rij|}}ij =Pj∈Nh(Ωi) ui−uj |rij |+ui−usym j |rsym ij | Pj∈Nh(Ωi)|ui−uj| |rij |+|ui−usym j| |rsym ij | (2.20) hPj∈Nh(Ωi)(ui−uj)+(ui−usym j) Pj∈Nh(Ωi)|ui−uj|+|ui−usym j|. Now, we use the following result for two sequences {ai}n i=1 {b}n i=1 (see [13] for further details): |Pn i=1 ai| Pn i=1 |ai|−|Pn i=1 bi| Pn i=1 |bi|=|Pn i=1 ai|−|Pn i=1 bi| Pn i=1 |ai|+ n X i=1 |bi|1 Pn i=1 |ai|−1 Pn i=1 |bi| ≤|Pn i=1 ai−bi| Pn i=1 |ai|+Pn i=1 |bi|−Pn i=1 |ai| Pn i=1 |ai|≤|Pn i=1 ai−bi|+Pn i=1 |ai−bi| Pn i=1 |ai| ≤2Pn i=1 |ai−bi| Pn i=1 |ai|.(2.21) Using simple algebraic manipulation, we have aq−bq= (a−b)Pq−1 k=0 akbq−ikfor q∈N+. For a, b ∈[0,1], it leads to |aq−bq| ≤ q|a−b|(see [13]). This inequality, together with
2.6. Lipschitz continuity 23 (2.20) and (2.21), leads to: 1 q|αi(u)−αi(v)|.Pj∈Nh(Ωi)((u−v)i−(u−v)j) + ((u−v)i−(u−v)sym j) Pj∈Nh(Ωi)|ui−uj|+|ui−usym j|.(2.22) On the other hand, the bounds |ui−uj| ≤ X k∈Nh(Ωi)|ui−uk|and |ui−uj| ≤ X k∈Nh(Ωj)|uj−uk|, (2.19), and (2.22), yield (νij(u)−νij(v))(ui−uj).qhd−1|v|X k∈Nsym h(Ωi)|(u−v)i−(u−v)k|(2.23) +qhd−1|v|X k∈Nsym h(Ωj)|(u−v)j−(u−v)k|. The second term is bounded by combining (2.18), (2.23), and the fact that the number of elements surrounding a node is bounded above independently of h: X i∈NhX j∈Nh(Ωi) (νij(u)−νij(v))`(i, j)ujzi.qhd−1|v||u−v|`|z|`. Next, we have to prove that the nonlinear mass matrix is also Lipschitz continuous. First, we note that X j∈Nh(Ωi) (1 −αi(uh))(ϕj, ϕi)uj+αi(uh)(1, ϕi)ui =X j∈Nh(Ωi) (ϕj, ϕi)uj+αi(uh)(ϕj, ϕi)(ui−uj). Thus hM(u)u, zi−hM(v)v, zi ≤ X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)(uj−vj)zi +X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)(ui−uj)(αi(uh)−αi(vh))zi +X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)((u+v)i−(u+v)j)αi(vh)zi.
24 Chapter 2. Monotonicity preserving stabilization for linear FEs Bounds for the second and third term follow the same lines as above. For the second term, we proceed as in (2.18), getting: X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)(ui−uj)(αi(uh)−αi(vh))zi .X i∈Nh 1 2X j∈Nh(Ωi) (ϕi, ϕj)(αi(uh)−αi(vh))2(ui−uj)2 1 2 ×kzk .qhd 2|u−v|`kzk. where we have used the spectral equivalence of the consistent and lumped mass matrices in the last inequality. The first and third term are easily bounded as X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)(uj−vj)zi≤ ku−vkkzk, X i∈NhX j∈Nh(Ωi) (ϕi, ϕj)((u+v)i−(u+v)j)αi(vh)zi≤qhd 2|u−v|`kzk. It proves the theorem. 2.7 Differentiable stabilization The previous nonlinear system is Lipschitz continuous, which improves the convergence of the nonlinear iterations. In fact, assuming that we supplement (2.1) with a diffusive term, existence and uniqueness can be proved in the diffusive regime (see [13]). However, even using Anderson acceleration nonlinear convergence can be very hard (see [64, 65] and Sect. 2.9). Based on these observations, we want to develop methods that lead to at least twice differentiable operators, i.e., ∂2T(uh) ∂2uh∈ C0, using the previous framework. This allows the usage of the Newton method to linearize the system, and reduces the required number of nonlinear iterations. Smoothness is achieved by substituting the non-differentiable functions of the previous formulation with smooth approximations. In order to end up with a twice differentiable method, we propose to use the following artificial diffusion: νij . = max σh{max σh{αεh,i(Fij(wh)), αεh,jFji(wh)},0},for i6=j, νii . =X j∈Nh(Ωi) j6=i νij.(2.24) The function max σh(·)is a regularized maximum function max σh{x, y}. =|x−y|1,σ 2+x+y 2,(2.25)
2.7. Differentiable stabilization 25 where |x|1,σ . =√x2+σis a smooth approximation of the absolute value. In order to keep dimensional consistency, σshould be a small parameter of order O|v|2`2(d−1), where `is a characteristic length of the problem. Let us define the smooth limiter function f(x)∈ C2that will be used in the definition of αεh, f(x). =(2x4−5x3+ 3x2+xif x < 1 1if x≥1. This function is used to smoothly limit the value of xup to 1. Further, let us define another smooth approximation of the absolute value, namely |x|2,ε . =x2 √x2+ε. Finally, the shock detector is defined as αεh,i(uh). = f Pj∈Nh(Ωi)J∇uhKij1,εh +γ Pj∈Nh(Ωi)2nn|∇uh·ˆ rij|2,εhooij +γ q ,(2.26) where γis a small parameter that prevents division by zero. It has been proved in Lemma 2.3.1 that αiequals 1 when iis an extremum in Ωi. Let us prove that this is still true for αεh,i. Lemma 2.7.1. If uhhas an extremum on i∈ Nhthen αεh,i(uh)=1. Proof. It is clear that f(x)equals 1 for x≥1, then the proof reduces to check that X j∈Nh(Ωi)J∇uhKij1,εh +γ≥X j∈Nh(Ωi) 2nn|∇uh·ˆ rij|2,εhooij +γ. Taking into account that px2+ε=|x|1,εh>|x| ≥ |x|2,εh=x2 √x2+ε,
26 Chapter 2. Monotonicity preserving stabilization for linear FEs and the fact that uj−uihas the same sign (or it is equal to zero) in all directions, it is easy to see that X j∈Nh(Ωi)J∇uhKij1,εh =X j∈Nh(Ωi) uj−ui |rij|+usym j−ui |rsym ij |1,εh ≥X j∈Nh(Ωi) uj−ui |rij|+usym j−ui |rsym ij | =X j∈Nh(Ωi) |uj−ui| |rij|+usym j−ui |rsym ij |≥X j∈Nh(Ωi) 2{{|∇uh·ˆ rij|}} ≥X j∈Nh(Ωi) 2nn|∇uh·ˆ rij|2,εhoo. It proves that αεh,i(uh) = 1 on an extremum. In fact, if the solution does not have an extremum, these quantities neither can have the same sign nor be zero in all cases. Since X j∈Nh(Ωi)J∇uhKij = lim ε→0X j∈Nh(Ωi)J∇uhKij1,εh and X j∈Nh(Ωi) 2{{|∇uh·ˆ rij|}} = lim ε→0X j∈Nh(Ωi) 2nn|∇uh·ˆ rij|2,εhoo, bound (2.8) leads to the fact that limε→0αεh,i(uh)<1when there is no extremum on i. It is straightforward to check the following results. Corollary 2.7.2. System (2.11) with the definition of the shock detector (2.26) and artificial diffusion (2.24) is LED and satisfies the local DMP. The method tends to a linearly preserving scheme as γ→0. Proof. From lemma 2.7.1 and the definition of the regularized maximum (2.25) it is easy to see that artificial diffusion in (2.24) is greater or equal to the one in (2.6). Hence, Theorem 2.4.2 still holds. The linearity preservation is straighforward. Remark 2.7.3. Note that the smoothed shock detector is not linearly preserving because αεh,i will never be zero. However, for regions where uhis constant the gradient is zero, thus the solution is not affected. In the case of uh∈P1(Ω), but not constant, αεh,i goes to zero with γ. Values of γof order 10−8(or even smaller) have been considered in the numerical experiments section with good nonlinear convergence properties. Thus, the linearity preservation is virtually preserved in practice. As in the previous section, when restricted to symmetric meshes, the following approximation (similar to the one in Barrenechea et al. [13]) of (2.26) maintains the same
2.8. Nonlinear Solvers 27 properties ˜αεh,i . = f Pj∈Nh(Ωi)ui−uj1,ε∗+γ∗ Pj∈Nh(Ωi)|ui−uj|2,ε∗+γ∗ q , with ε∗∼ O(h2ε)and γ∗∼ O(hγ). 2.8 Nonlinear Solvers In this section the methods used for solving the system of nonlinear equations resulting from the above formulation (2.12) with the artificial diffusion defined in (2.24) is discussed. Taking advantage of the differentiability of the stabilization described in Sect. 2.7, Newton’s method is used for the smooth version of the method. In addition, we use fixed point iterations with Anderson acceleration to compare against Newton’s method performance. In order to define the schemes, it is useful to write the time-discrete problem (2.12) as A(un+1 h)un+1 h=G where Gis the force vector. Let J(un+1 h). =∂T(un+1 h) ∂un+1 h be the Jacobian. Since the above problem is nonlinear we will solve it iteratively. We denote by uk,n+1 h the k-th iteration of uhat time step n+ 1. Let us define some auxiliary variables used in the definition of the algorithms: mdenotes the number of previous nonlinear iterations used in Anderson acceleration, sis the slope resulting form fitting the last m nonlinear errors, smin is the minimum slope allowed before increasing the relaxation, ωis the relaxation parameter, ωmin is its allowed minimum, kmax is the maximum nonlinear iterations allowed, tol is the nonlinear tolerance, and nlerr is the nonlinear error. For the non-differentiable methods in Sect. 2.3 we use Picard linearization with Anderson acceleration (see Alg. 1). Our particular implementation also includes a simple convergence rate test, where it is decided if the relaxation parameter should be reduced or not. This improves the global convergence rate and the robustness of the method. Moreover, we add a projection onto Vadm hto ensure that the global DMP in Def. 2.5.1 is satisfied at all nonlinear iterations. This step is of special interest in the case of solving for variables that cannot become negative, e.g., the density. In this case, the projection onto the space of admissible solutions is performed truncating the obtained solution. However, more sophisticated methodologies can be also applied but at a higher computational cost. For the differentiable method, Newton’s linearization is used (see Alg. 2). In addition, we supplement it with the line search method to improve robustness. We use numerical 1D minimization of the residual norm up to a tolerance of 10−4for the line search method. Following the same approach in Alg. 1, a projection to the FE space of admissible solutions can be performed in Alg. 2. As said before, this step ensures that for all nonlinear iterations the solution satisfies the global DMP. The numerical experiments
28 Chapter 2. Monotonicity preserving stabilization for linear FEs in the next section show that the modified method keeps quadratic convergence, even though we do not have a theoretical analysis. Algorithm 1: Fixed point iterations with relaxed Anderson acceleration Input:u0,n+1 h,m,smin,ωmin,tol,A,G,kmax Output:uk,n+1 h,k k= 1,nlerr1=tol while (nlerrk≥tol) and (k < kmax)do Set mk= min(k, m) Solve A(uk,n+1 h)˜uk,n+1 h=G Compute rk,n+1 = ˜uk,n+1 h−uk,n+1 h Minimize kPmk i=1 ξk irk−mk+i,n+1kwith respect to ξk isubject to Pmk i=1 ξk i= 1 Set uk+1,n+1 h= (1 −ωk)Pmk i=1 ξk iuk−m+i,n+1 h+ωkPmk i=1 ξk i˜uk−mk+i,n+1 h Project uk+1,n+1 hto Vadm h Set nlerrk=kuk+1,n+1 h−uk,n+1 hk kuk+1,n+1 hk Compute the slope (s) of {nlerri}with k≥i≥k−mk if (s<smin)and (ω > ωmin)then Set ωk+1 =ωk−0.1 else Set ωk+1 =ωk Update k=k+ 1 Algorithm 2: Newton’s method + Line search Input:u0,n+1 h,un h,tol,J,R,kmax Output:uk,n+1 h,k k= 1,nlerr1=tol while (nlerrk≥tol) and (k < kmax)do Solve J(uk,n+1 h)∆uk,n+1 h=−T(uk,n+1 h) Minimize kT(uk,n+1 h+ξk∆uk,n+1 h)kwith respect to ξ∈[0,1] Set uk+1,n+1 h=uk,n+1 h+ξk∆uk,n+1 h Project uk+1,n+1 hto Vadm h Set nlerrk=kξk∆uk,n+1 hk kuk+1,n+1 hk Update k=k+ 1 2.9 Numerical Experiments 2.9.1 Steady problems First, in order to test the previous formulation, the convergence to a smooth solution is analyzed. For this purpose, the following equation is solved ∇·(vu) = 0 in Ω = [0,1] ×[0,1], u=uDon Γin,(2.27)
2.9. Numerical Experiments 29 with v(x, y). = (1,0), and inflow boundary conditions uD=y−y2on ∂Ω\{x= 1}. This problem consists in the transport of the parabolic profile along the xdirection, which has the analytical solution u(x, y) = y−y2. Fig. 2.2 shows the convergence rates using the previously defined formulation ((2.12) with (2.24)), and the Galerkin formulation. To perform this test, an initial mesh of 12 ×12 Q1has been considered, then successive refinements have been performed up to a96 ×96 Q1mesh. Analogous meshes has been also used for P1FE. Newton’s method has been used with q= 4,ε= 10−7,σ=|v|h410−8and γ= 10−10. In this case, σhas been scaled as |v|2L2(d−3)h4in order to recover optimal convergence, where Ldenotes a characteristic length of the physical domain Ω. As desired, the convergence rates are not affected by the stabilization, while (as expected) the stabilized solutions have higher errors. nx=ny 20 30 40 50 60 70 80 90 kuh!uk2 10-5 10-4 10-3 10-2 P1 Galerkin Q1 Galerkin P1 Stabilized Q1 Stabilized Figure 2.2: Convergence test, L2(Ω) error versus size of the mesh. For P1and Q1FE meshes ranging from h= 1/12 to h= 1/96. Newton’s method has been used with parameters q= 4,ε= 10−7,σ=|v|h410−8 and γ= 10−10. A typical linear test to assess the performance of a shock capturing method is the propagation of a discontinuity. Consider now the previous hyperbolic PDE (2.27) with v(x, y). = (1 /2,sin −π/3), and inflow boundary conditions uD= 1 on {x= 0}∩{y > 0.7} and y= 1, while uD= 0 at the rest of the inflow boundary. This problem has the following analytical solution u(x, y) = (1if y > 0.7+2xsin −π/3, 0otherwise. At Fig. 2.3(a), the numerical solution using the stabilization in (2.24) is shown. A 48 ×48 Q1mesh have been used. The values chosen for the parameters in (2.24) are q= 25,ε= 10−4,σ=|v|10−9, and γ= 10−10. This parameter choice makes the solution at the outflow sharp while the DMP is always satisfied. Furthermore, convergence is not
30 Chapter 2. Monotonicity preserving stabilization for linear FEs jeopardized thanks to the smoothed stabilization. Particularly, it took 18 iterations for the Newton’s method to converge to a nonlinear tolerance of 10−6. The non-smooth version in Fig. 2.3(b) ((2.11) with (2.6)) did not converge using Anderson acceleration, adding a fixed relaxation parameter of ω= 0.5took 392 iterations, and 117 with Alg. 1. In any case, observing Fig. 2.4, where the outflow profile is depicted, no apparent improvement on accuracy is observed when using the non-smooth version. (a) Smooth stabilization (2.26), with parameters q= 25,ε= 10−4,σ= |v|10−9, and γ= 10−10. (b) Non-smooth version (2.7) with q= 25. Figure 2.3: Stabilized solution of the straight propagation of a discontinuity test using the steady version of discrete problem (2.12) with two stabilization choices (2.26) or (2.7). x 0 0.2 0.4 0.6 0.8 1 u 0 0.2 0.4 0.6 0.8 1 Non-smooth version Smoothed stabilization Figure 2.4: Stabilized solution of the straight propagation of a discontinuity test using the steady version of discrete problem (2.12) with two stabilization choices (2.26) and (2.7). The stabilization parameters used for the smoothed version are q= 25,ε= 10−4,σ=|v|10−9, and γ= 10−10. Fig. 2.5 shows the solution for several combinations of qand ε, with σ=|v|ε10−5 and γ= 10−10, solved with the two nonlinear solvers presented in the previous section over a 48 ×48 Q1mesh. Furthermore, the ku−uhkL1and ku−uhkerrors, computed at the whole domain and restricted to the outflow boundary, are listed in Table 2.1. These results show that, as expected, either increasing qor reducing εthe L1error diminishes. Nevertheless, the computational cost also increases at a higher rate. The same can be observed for the 2error. It is slightly reduced after increasing qor diminishing ε, while this makes nonlinear convergence much harder. Moreover, comparing both nonlinear solvers in Sect. 2.8, it is important to note that using Newton’s method the number of nonlinear iterations is reduced between 10 to 15 times.
2.9. Numerical Experiments 37 (a) Initial conditions. (b) LED scheme. (c) Global DMP scheme. (d) LED DMP nonsmooth stabilization. Figure 2.12: 3 Body rotation test results using discrete problem (2.12) and two different artificial diffusions ((2.24) leading an LED scheme, and (2.15) with (2.26) leading a global DMP scheme). Using a 150 ×150 Q1 element mesh, and parameters: q= 25,γ= 10−8,σ=|v|10−10,ε= 10−4, and ∆t= 10−3. 2.9.3 Burgers’ equation Finally, let us test our stabilization with a nonlinear transient problem. Particularly the 2D Burgers’ equation, i.e. equation (2.28) with v. = (1,1)u /2, is solved on Ω = [0,1]×[0,1] using a 150×150 Q1mesh. The discretization in time is performed using the BE method
38 Chapter 2. Monotonicity preserving stabilization for linear FEs with a time step of 10−2. The initial conditions at t= 0 are u0(x, y). = −0.2if x < 0.5and y > 0.5 −1if x > 0.5and y > 0.5 0.5if x < 0.5and y < 0.5 0.8if x > 0.5and y < 0.5 , and the solution is advanced until t= 0.5. (a) Solution for:q= 1,ε= 10−2,σ= |v|10−6, and γ= 10−8. (b) Solution for: q= 4,ε= 10−3,σ= |v|10−7, and γ= 10−8. Figure 2.13: Burger’s equation solutions at t= 0.5using discrete problem (2.12) and (2.6) with (2.24). Using a 150 ×150 Q1element mesh, ∆t= 10−2, and two sets of parameters q,γ,σ, and ε. The following stabilization parameters have been used for obtaining the results in Fig. 2.13(a): q= 1,ε= 10−3,σ=|v|10−6, and γ= 10−8. Although the parameters used are not enforcing a particularly sharp solution (see Figs. 2.5 and 2.8), Fig. 2.13(a) shows properly transported and minimally smeared shocks. Only in the lower right region the method appears to be more diffusive than desired. Notice that in that region the gradient in the xdirection spreads as yincreases, while it should not. Nevertheless, in Fig. 2.13(b), that shows the solution for q= 4,ε= 10−4,σ=|v|10−7, and γ= 10−8. the method is less diffusive and the obtained shocks are even sharper. In any case, both choices satisfy the DMP for all time steps. 2.10 Conclusions In this chapter, we have considered a nonlinear stabilization technique for the FE approximation of scalar conservation laws with implicit time stepping. The method relies on an artificial diffusion method, based on a graph-Laplacian operator. The artificial diffusion is judiciously chosen in order to satisfy a local DMP for steady problems. It is nonlinear, since it depends on a shock detector. Further, the resulting method is linearity preserving. The same shock detector is used to gradually lump the mass matrix.
2.10. Conclusions 39 The resulting method is LED, positivity preserving, and also satisfies a global DMP. Lipschitz continuity has also been proved. However, the resulting scheme is highly nonlinear, leading to very poor nonlinear convergence rates, even using Anderson acceleration techniques. It is due to the fact that the nonlinear operator to be inverted at every time step is non-differentiable. The critical problem of nonlinear convergence of implicit monotonic methods based on nonlinear artificial diffusion have already been previously reported in the literature (see [57]). As a result, we propose a smooth version of the scheme. It leads to twice differentiable nonlinear stabilization schemes, which allows one to straightforwardly use Newton’s method using the exact Jacobian. Twice differentiability ensures quadratic convergence. We have considered two nonlinear solvers, namely Anderson acceleration and Newton’s method. We have observed numerically that the effect of the smoothness has a positive impact in the reduction of the computational cost. The impact of using Newton’s method versus Anderson acceleration is also very positive. In general, using the Newton method with a smooth version of the method we can reduce 10 to 20 times the number of iterations of Anderson acceleration with the original non-smooth algorithms. All the monotonic properties are satisfied (as theoretically proved) in the numerical experiments. Steady and transient linear transport, and transient Burgers’ equation have been considered in 2D. In any case, these properties are only true for the converged solution, but not for iterates. In this sense, we have also proposed the concept of projected nonlinear solvers, where a projection step is performed at the end of every nonlinear iterations onto a FE space of admissible solutions. The space of admissible solutions is the one that satisfies the desired monotonic properties (maximum principle or positivity). The projection has no effect on the quality of the nonlinear convergence. Future work should tackle the entropy stability analysis of the resulting schemes when applied to nonlinear problems. Some initial results in this direction can be found in [23]. The extension to systems of conservation laws and higher order methods in space and time is another interesting line of research.
Chapter 3 Arbitrary order space–time monotonicity preserving scheme This chapter is devoted to a nonlinear stabilization technique for convection–diffusion– reaction and pure transport problems discretized with space–time isogeometric analysis. The stabilization is based on a graph-theoretic artificial diffusion operator and a novel shock detector for isogeometric analysis. Stabilization in time and space directions are performed similarly, which allow us to use high-order discretizations in time without any CFL-like condition. The method is proved to yield solutions that satisfy the discrete maximum principle (DMP) unconditionally for arbitrary order. In addition, the stabilization is linearity preserving in a space–time sense. Moreover, the scheme is proved to be Lipschitz continuous ensuring that the nonlinear problem is well-posed. Solving large problems using a space–time discretization can become highly costly. Therefore, we also propose a partitioned space–time scheme that allow us to select the length of every time slab, and solve sequentially for every subdomain. As a result, the computational cost is reduced while the stability and convergence properties of the scheme remain unaltered. In addition, we propose a twice differentiable version of the stabilization scheme, which enjoys the same stability properties while the nonlinear convergence is significantly improved. Finally, the proposed schemes are assessed with numerical experiments. In particular, we considered steady and transient pure convection and convection–diffusion problems in one and two dimensions. 3.1 Introduction Many different applications in science and industry require solving problems satisfying some sort of positivity or maximum principle (MP) property. These include scalar transport problems, compressible flows, or fluid-based MHD simulations, among others. These problems are of particular interest in a variety of industries and scientific research areas, such as the chemical industry, aviation, aerospace, or nuclear fusion research, just to cite few examples. Some of these problems exhibit a multiscale nature in time. In those cases, explicit methods are not suitable, since the smallest time scales pose very stringent stability conditions to the time step length, i.e., fully resolved time simulations are required. 41
42 Chapter 3. Arbitrary order space–time monotonicity preserving scheme Thus, implicit methods are favored in applications where the smallest time scales are not of scientific or engineering interest. As a result, schemes that preserve monotonicity (or at least positivity) for implicit time integration are of special interest. The standard technique to attain such schemes is adding nonlinear artificial diffusion (usually called shock capturing). The common ingredients of a shock capturing or nonlinear stabilization method are the following. The first ingredient is the artificial diffusion, which needs to be sufficient to eliminate non-physical oscillations. The schemes in [7, 8, 23, 25] use an element-based artificial diffusion with a standard PDE-based diffusion operator. The drawback of this choice is the fact that the DMP only holds under unpractical mesh restrictions. This problem has been solved by Guermond and Nazarov in [38, 41] by replacing the PDE-based diffusion operator by an edge or graph-theoretic diffusion operator; see [4, 5, 65, 71, 75] for schemes that preserve the DMP on arbitrary meshes using a graph-Laplacian. The second ingredient is a shock detector, which is the term responsible of deactivating the artificial diffusion in smooth regions. A good shock detector is of vital importance for minimizing the numerical diffusion while satisfying a DMP. One example of shock detector is the one developed in 1D by Burman in [23] and later extended to multiple dimensions by Badia and Hierro [7]. The last ingredient consists on perturbing the mass matrix. One option is a full lumping of the mass matrix, but it can lead to unacceptable phase errors. Instead, a nonlinear lumping is used, e.g., in [4, 5], using the same shock detector to lump the mass matrix. Other alternatives can be found in [40, 65]. It is worth mentioning that all previous stabilization methods yield a very stiff nonlinear system of equations. In fact, some of the methods proposed in the literature are not even Lipschitz continuous and thus ill-posed (see [13]). In practice, the nonlinear convergence of these methods is unacceptably slow, making hard its practical use. To solve this problem, Badia et al. [4, 5] have designed differentiable nonlinear stabilization terms, noticeably improving the nonlinear convergence. The methods commented above have an algebraic nature and provide some type of DMP for the nodal values. The monotonicity of the nodal values only translates into monotonic solutions if the FE space satisfies the convex hull property, which is only true in the first order case. As a result, using the ideas above it does not seem possible to design monotonic second or higher order methods. Recently, Kuzmin and coworkers [2, 71], have proposed instead the usage of Bernstein–Bèzier FEs, since they satisfy the convex hull for high-order. However, the temporal dimension is discretized using Backward Euler or SSP RK methods (see [53]). In the first case, the problem is first order in time, whereas in the second case, a CFL-like condition arises [67], since high-order SSP methods pose a restriction on the time step size similar to the ones in explicit methods [53]. The main contribution present in this chapter is the development of a high-order (both in space and time) and DMP-preserving discretization for the convection–diffusion– reaction and pure transport problems. This is achieved by combining the nonlinear
3.2. Preliminaries 43 stabilization techniques in the previous chapter and [5] with a new shock detector for arbitrary order space–time isogeometric analysis. Another novelty introduced in this chapter is the stabilization in the time direction, which is performed in a similar manner as in space. This results in an unconditionally stable high-order method in time (and space). However, the space–time method requires to solve the whole space–time problem at once, which increases the computational cost. Hence, we also propose a partitioned approach in the temporal direction, where one can determine the width of the time slab to be computed every time. This strategy allows us to maintain a reasonable computational cost while having a high-order scheme in space and time, as well as satisfying the DMP without any CFL-like condition. Finally, we also propose a differentiable version of the above scheme. This allow us to use Newton’s method, which improves nonlinear convergence significantly. This chapter is structured as follows. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems in Sect. 3.2. Then, the stabilization techniques are introduced in Sect. 3.3. Sect. 3.4 is devoted to the partitioned time integration scheme. Afterwards, we introduce a regularized version of the stabilization term in Sect. 3.5. Finally, we show numerical experiments in Sect. 3.6 and draw some concluding remarks in Sect. 3.7. 3.2 Preliminaries 3.2.1 Convection–Diffusion problem We consider a transient convection–diffusion problem with Dirichlet boundary conditions. Let Ω×(0, T). =Qd+1 α=1(0, Lα)be a (d+ 1)-cube, where dis the number of spatial dimensions. Then, the problem reads: ∂tu+∇·(vu)−∇·(µ∇u) = gin Ω×(0, T], u(x, t) = u(x, t)on ∂Ω×(0, T], u(x, 0) = u0(x)x∈Ω, (3.1) where vis a divergence-free convection velocity, µ≥0is a scalar constant diffusion, and g(x, t)is the body force. In the case of pure convection (µ= 0), boundary conditions are only imposed at the inflow Γin . ={x∈∂Ω : v·n∂Ω<0}, where n∂Ωis a unit vector outward-pointing normal to the boundary. We also define the outflow boundary as Γout . =∂Ω\Γin. Moreover, we will also consider the steady problem, which is obtained by dropping the time derivative term and the initial condition. It is important to mention that a reaction term can be included without harming any of the properties satisfied by the schemes introduced below. However, a convection–diffusion–reaction problem only satisfies a MP if the minimum is negative and the maximum positive (analogously for its proposed discretizations). In other words, it only satisfy a weak MP, see [13]. In order
44 Chapter 3. Arbitrary order space–time monotonicity preserving scheme to simplify the discussion below, we will limit the present chapter to pure convection and convection–diffusion problems. In order to avoid technicalities and facilitate the exposition of the stabilization method, we restrict this chapter to cubic domains. However, it is possible to work with complex geometries using standard procedures from isogeometric analysis [29]. E.g., a complex geometry would be divided in several parts, which would be mapped to multiple d-cubic patches. The stabilization method presented in this chapter is independent from this procedure. 3.2.2 Discretization In this chapter, we consider a standard B-spline discretization with interpolative boundaries (see [29]). A spline of order pin the variable xis a piecewise polynomial function in xof degree p. The values of xin which different polynomials meet are call knots. Knots might be placed at the same location, i.e. can be repeated. When the knots are not repeated, the first p−1derivatives of the spline are continuous. When a knot is repeated rtimes, only the first p−rderivatives are continuous across that knot. Knots are sorted in increasing order and collected in the so called knot vector {ξ1, ξ2, . . .}. Given a knot vector, B-splines of order pare basis functions for spline functions of the same order. B-splines are constructed in a recursive way using the Cox-de Boor formula: B0 i(x). =(1if ξi≤ξ < ξi+1 0otherwise , Bk i(x). =x−ξi ξi+k−ξi Bk−1 i(x) + ξi+k+1 −x ξi+k+1 −ξi+1 Bk−1 i+1 (x), for k= 1, ..., p. By construction, Bp i(x)has compact support, is non-negative, and nonzero in [ξi, ξi+p+1]. Notice that its support increases with the degree of the polynomial. Let us consider the domain [0, L]and the uniform partition into msub-intervals of size h=L/m. The open knot vector {ξ1, . . . , ξm+2p+1}is defined as follows. The first p+ 1 knots are located at zero, i.e., ξ1=. . . =ξp+1 = 0. The last p+ 1 knots are located at L, i.e., ξm+p+1 =. . . =ξm+2p+1 =L. The interior points are equidistributed, with ξi= (i−p−1)h, for i=p+ 1, . . . , m +p+ 1. It leads to a basis Bp i(x)(for i= 1, . . . , m +p) for a space of splines in [0, L]and a partition of unity, i.e., Pm i=1 Bp i(x)=1for x∈[0, L]. Any spline v(x)of order pin [0, L]can uniquely be defined by the control points (v1, . . . , vm+p)∈Rm+das the linear combination of B-splines v(x) = Pm+p i=1 Bp i(x)vi. In one dimension, the basis functions obtained from an open knot vector are interpolatory at the extremes, i.e., v(0) = v1and v(L) = vm+p+1 (see Fig. 3.1). For a first order polynomial vin [0, L], it holds v(x) = Pm i=1 Bp i(x)v(xi), where xi. = (ξi+1 +. . . +ξi+p)/p are called the Greville abscissae [30, 76]. Let us consider the number of partitions per dimension with mα, for α= 1, . . . , d + 1. We represent with Nhthe set of multi-indices i. = (i1, ..., id+1)∈Zd+1 with iα∈
3.2. Preliminaries 45 {1, . . . , mα+p}. Every i∈ Nhcan be expressed as (ix, it), where ixis the spatial index and itis the temporal index. The (d+ 1)-dimensional B-spline is defined as the tensor product of d+ 1 unidimensional B-splines Bp i(x). =Bp i1(x1)×···×Bp id+1 (xd+1). Notice that a Greville abscissa in the case of a multidimensional spline reads xi= (xi1, ..., xid+1 ). We define the space of splines Vh. = span{Bp i(x) : i∈ Nh}. We use the notation ϕi≡Bp i. The order is omitted since it is assumed to be fixed. Thus, every spline vh∈Vh can be written as vh=Pi∈Nhϕivi. Furthermore, we define the following sets of indices, which are useful for the definition of the forthcoming schemes. The set of neighbors of i is defined as Ni h . ={j∈ Nh:|i−j|∞≤1}. We define as Si h . ={j∈ Nh:|i−j|∞≤p} the set of indices whose associated shape functions intersect with the support of ϕi. Figure 3.1: Representation of the basis functions of V2 hin one dimension, with its associated Greville abscissae. We use standard notation for Sobolev spaces. The L2(ω)scalar product is denoted by (·,·)ωfor ω⊂Ω. However, we omit the subscript for ω≡Ω. The L2(Ω) norm is denoted by k·k. 3.2.3 Discrete problem The weak form of (3.1) using the Galerkin method reads: find uh∈Vhsuch that uh(x, t) = uh(x, t)on ∂Ω×(0, T],uh(x,0) = u0h(x)on Ω×{0}, and (∂tuh, vh)+(v·∇uh, vh) + µ(∇uh,∇vh)=(g, vh),∀vh∈Vh,(3.2) where uh(t)and u0hare projections of u(t)and u0to Vh, respectively, such that the local DMP is satisfied (see Def. 3.2.2). Furthermore, we can rewrite the previous discrete problem in matrix form as Kijuj=Fi, where Kij . = (∂tϕj, ϕi)+(v·∇ϕj, ϕi)+µ(∇ϕj,∇ϕi), and Fi. = (g, ϕi)for i,j∈ Nh. Notice that we have not applied the boundary conditions yet. To apply boundary conditions the space of test functions is restricted to vh∈Vh0, and the force vector is redefined as Fi. = (g, ϕi)−(∂tuh, ϕi)−(v·∇uh, ϕi)−µ(∇uh,∇ϕi). 3.2.4 Monotonicity properties In this section we define all the properties that we demand our scheme to fulfill. In this case, since we are using a space–time discretization, it becomes more useful to define
46 Chapter 3. Arbitrary order space–time monotonicity preserving scheme these properties in a space–time sense. This means that the variation of uhin the temporal direction will also be taken into account to define an extremum. Hence, we define the concept of a local discrete extremum as follows. Definition 3.2.1 (Local Discrete Extremum).The function uh∈Vhhas a local discrete minimum (resp. maximum) on i∈ Nh, if ui≤uj(resp. ui≥uj)∀j∈ Ni h. For problems that satisfy a maximum principle, e.g., problem (3.2) with g= 0, it is also important to define the concepts of local and global space–time DMP. The latter is a slightly weaker property than the former, but it is more useful for the discussion in this chapter. A local DMP is a stronger property because it implies that no oscillations can appear, while the global DMP only implies that the global extrema are located at the boundary conditions. Definition 3.2.2 (Local space–time DMP).A solution uh∈Vhsatisfies the local discrete maximum principle if for every i∈ Nh min j∈Ni h\{i}uj≤ui≤max j∈Ni h\{i}uj. Definition 3.2.3 (Global space–time DMP).A solution uh∈Vhsatisfies the global discrete maximum principle if the global extrema are located at boundary conditions, i.e., for every i∈ Nh min min x∈∂Ω, t∈[0,T) uh(x, t),min x∈Ωuh0(x) ≤ui≤max max x∈∂Ω, t∈[0,T) uh(x, t),max x∈Ωuh0(x) . Finally, let us recall the definition of linearity-preservation, which is a desired property to achieve high-order convergence in smooth regions (see [66]). Definition 3.2.4 (Linearity-preservation).A stabilization term, Bij(uh), is said to be linearity-preserving if, for a solution that is linear in all directions in the neighborhood of xi, then the stabilization term becomes null, i.e., Bij(uh) = 0 if uh(x)∈ P1(Ωi)where Ωiis the convex hull defined by the set of neighboring Greville abscissae {xj}j∈Nk h,k∈Si h. 3.3 Lipschitz-continuous nonlinear stabilization In this section we define a nonlinear stabilization operator, Bh(wh;uh, vh), to be added to the discrete problem (3.2), such that it satisfies at least the global DMP in Def. 3.2.3. Let us define Bij(uh). =Bh(uh;ϕj, ϕi). We also enforce that, for any uh∈Vh,Bij(uh) 1. has compact support: Bij (uh) = 0 if j6∈ Si h, 2. is symmetric: Bij(uh) = Bji(uh), 3. is conservative: Pj∈Si h\{i}Bij(uh) = −Bii(uh).
3.6. Numerical experiments 53 3.6 Numerical experiments In this section we present some numerical experiments showing the behavior of the scheme previously introduced. First, a convergence analysis is performed in order to assess the correctness of the proposed scheme and its implementation. Then, we assess the performance of the proposed stabilization method for high-order discretizations, including a brief analysis of the effect of the regularization. 3.6.1 1D Transient Diffusion The purpose of this test is assessing the partitioned time integration scheme in Sect. 3.4. To this end, we solve the following problem for t∈(0,1] and x∈Ω. = (0,1), (∂tu+∂xxu=fin Ω×(0,1] u= 0 at ∂Ω,(3.11) where f. = 2(6x2−6x+ 1)(t(t−1))2+ 2t(t−1)(2t−1)(x(x−1))2. This problem has u= (x(x−1))2(t(t−1))2as exact solution. We perform a convergence analysis where the mesh is successively refined in the time direction for first, second, and third order discretizations. In particular, the distance between knots in the temporal direction is refined as δt ={0.2,0.1,0.05,0.025,0.0125}for the first order discretization. In spatial directions, the distance is small enough (h= 1/400) to prevent that spatial discretization errors affect the analysis. Second and third order discretizations are obtained using the following k-refinement (see [29] for more details). We refine the discretization such that the number of control points increase at the same rate as a Lagrangian FE discretization does when its order is increased. Fig. 3.3 shows the result of k-refinements to p= 2 and p= 3 discretizations, for an interior subset of the discretization. Henceforth, we will use this kind of k-refinement in order to increase the discretization order. Figure 3.3: Second and third order discretizations obtained from the k-refinement of an initial first order discretization. Notice that shape functions are depicted for interior knots, at boundary knots shape functions become interpolatory, see Fig. 3.1. We measure the relative L2norm and H1semi-norm of error in the whole space– time domain, and compute the resulting convergence rate. Errors in L2norm and H1 semi-norm are depicted in Fig. 3.4 (a) and (b), respectively. In Table 3.1 the measured convergence rates are shown for the original non-partitioned scheme and the proposed in
54 Chapter 3. Arbitrary order space–time monotonicity preserving scheme 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 log10 1/∆t −12 −10 −8 −6 −4 −2 log10 ku−uhkL2/kuhkL2 1 2 1 3 1 4 1st order p-ST 1st order ST 2nd order p-ST 2nd order ST 3rd order p-ST 3rd order ST (a) Relative L2norm of the error. 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 log10 1/∆t −10 −8 −6 −4 −2 0 log10 |u−uh|H1/|uh|H1 1 1 1 2 1 3 1st order p-ST 1st order ST 2nd order p-ST 2nd order ST 3rd order p-ST 3rd order ST (b) Relative H1norm of the error. Figure 3.4: Convergence in time results for problem (3.11), using standard and partitioned space-time schemes. Sect. 3.4. We observe a slight increase in the error for the partitioned scheme. However, the obtained results show optimal convergence rates for both schemes introduced above. Table 3.1: Measured convergence rates in L2norm and H1semi-norm, for problem (3.11). Order Method L2convergence H1convergence 1 p-ST -1.98 -0.98 1 ST -2.04 -0.96 2 p-ST -3.03 -2.00 2 ST -3.00 -2.01 3 p-ST -3.99 -3.00 3 ST -3.99 -2.99 p-ST: Partitioned space-time, ST: space-time. 3.6.2 Steady convection In this experiment we assess the convergence of the stabilized schemes introduced in Sect. 3.5. We use a steady pure convection problem with a non-monotonic smooth solution. In particular, we solve the following problem for x∈Ω. = [0,1]2, (v·∇u= 0 in Ω u=uDat Γin ,(3.12) where uD= sin 2πx−y tan θ,v= (cos θ, sin θ), and θ=π/3. The analytical solution of the above problem reads u= sin 2πx−y tan π/3. The convergence analysis is performed for first, second, and third order discretizations, i.e. p={1,2,3}. We use a standard nonlinear solver (see [5] for details), with a nonlinear tolerance uk+1−uk uk<10−6. The following stabilization parameters have been used: q= 10,ε= 10−5,σ= 10−6, γ= 10−10. The selection of these parameters is based on the outcome of previous works [4, 5] and Sect. 3.6.3.
3.6. Numerical experiments 55 Convergence plots are shown in Fig. 3.5 and the corresponding convergence rates in Table 3.2. As expected, it is observed that the scheme recovers second order convergence in the L2error norm and first in the H1error semi-norm. It is known that the stabilized scheme should recover second order convergence for p= 1. However, due to peak clipping errors, higher convergence rates are not expected even if a higher order discretization is used [58]. In any case, we do observe that the error diminishes as the discretization order is increased using the k-refinement previously defined. 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 log10 1/∆x −3.5 −3.0 −2.5 −2.0 −1.5 log10 ku−uhkL2/kuhkL2 1 2 1st order 2nd order 3rd order (a) Relative L2norm of the error. 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 log10 1/∆x −1.0 −0.8 −0.6 −0.4 −0.2 0.0 0.2 log10 |u−uh|H1/|uh|H1 1 1 1st order 2nd order 3rd order (b) Relative H1semi-norm of the error. Figure 3.5: Convergence in space results for problem (3.12). Table 3.2: Measured convergence rates in L2and H1norms, for problem (3.12). Order L2convergence H1convergence 1 1.77 0.88 2 1.85 0.96 3 1.86 0.99 3.6.3 Nonlinear convergence In the current test, we aim to briefly analyze the effect of the stabilization parameters on the nonlinear convergence of the method. To this end, we solve the following 1D pure convection problem with discontinuous initial conditions. ∂tu+v·∇u= 0 in Ω×[0, T) u=u0at t= 0 u=uDat ∂Ω ,(3.13) where v. = 1,Ω. = (0,1],T= 0.5, and u0. = 1−H0(x−0.25), where H0is the well-known zero-centered Heaviside function. First and second discretization orders are used in a coarse mesh of 25 ×25 control points. To obtain the second order mesh, we perform the k-refinement as in the previous experiment.
56 Chapter 3. Arbitrary order space–time monotonicity preserving scheme We refer the reader to [4, 5] for a deeper analysis on the effect of each regularization parameter. Therein, the same family of shock detectors is used in the context of first order cG and dG Lagrangian FEs. In this chapter, we analyze the effect of the regularization globally using a fixed relation between the different parameters. In particular, we use the following parameters: γ= 10−10,σ=ζ,ε=ζ2, where ζ={10−1,10−2,10−3,10−4}. Furthermore, the effect is also compared as qis incremented, particularly for q={1,2,5,10}. In addition, the non-regularized version is also used to show the improvement in the nonlinear convergence. The relaxed Picard and hybrid nonlinear solvers presented in [5] are used, and the nonlinear tolerance is set to 10−5. 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 34; 17 (38) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 36; 17 (38) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 37; 19 (38) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 36; 19 (38) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 45; 20 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 53; 18 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 50; 23 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 50; 23 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 46; 33 (81) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 63; 57 (81) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 79; 57 (81) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 72; 64 (81) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 49; 38 (136) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 100; 57 (136) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 97; 51 (136) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 85; 59 (136) ζ= 10−1ζ= 10−2ζ= 10−3ζ= 10−4 q= 10 q= 5 q= 2 q= 1 Figure 3.6: Effect of the regularization parameters for first order discretizations. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picard scheme. Fig. 3.6 and 3.7 show the results for first and second order discretizations, respectively. In general terms, as qis increased or ζis decreased, sharper solutions are observed. However, nonlinear iterations increase. As expected, the hybrid method outperforms the relaxed Picard method. Even though it requires more nonlinear iterations, the nonregularized detector might be a simpler (parameter-free) alternative to the regularized one. Finally, it is worth mentioning that a slight increase in the required number of iterations is observed as the discretization order is increased. However, the obtained results are more accurate, i.e., the discontinuity becomes sharper.
3.6. Numerical experiments 57 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 52; 22 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 48; 22 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 54; 23 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 48; 26 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 52; 24 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 52; 31 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 51; 32 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 52; 31 (53) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 63; 32 (67) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 56; 28 (67) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 58; 41 (67) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 62; 43 (67) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 81; 25 (91) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 111; 41 (91) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 98; 36 (91) 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 75; 53 (91) ζ= 10−1ζ= 10−2ζ= 10−3ζ= 10−4 q= 10 q= 5 q= 2 q= 1 Figure 3.7: Effect of the regularization parameters for second order discretization. The numbers in legends are the number of nonlinear iterations performed. First number is for relaxed Picard and the next for hybrid scheme, both for the regularized stabilization. The number in brackets is the number of iterations required to converge the non-differentiable method using relaxed Picard scheme. 3.6.4 1D Sharp layer propagation The performance of the stabilization schemes is analyzed as the discretization order is increased. To this end, we use again the previous problem (3.13). The regularization parameters are kept fixed, while the discretization is modified both in terms of the order of accuracy and the number of control points. We use a nonlinear tolerance of 10−5. The regularization parameters used are q= 10, ε= 10−8,σ= 10−6, and γ= 10−10. With this setting, we solve the above problem using a discretization that keeps the number of control points fixed as the order is increased, and another one using the k-refinement defined in the previous experiment. For the former, we use a discretization of 120 by 60 control points. For the latter, we start with a first order discretization of 120 by 60 control points and refine as previously explained. In Fig. 3.8(a) the solution at t= 0.5is shown for different orders and fixed number of control points, and using the k-refinement in Fig. 3.8(b). We observe that for nonsmooth solutions, fixing the number of control points and increasing the order does not improve the results. This is a consequence of the underlying discretization properties. The support of the shape functions becomes larger as the order is increased. Therefore, nonsmooth solutions become slightly more smeared. In the case of Fig. 3.8(b), as expected, we observe better approximations as the order is increased using the k-refinement.
58 Chapter 3. Arbitrary order space–time monotonicity preserving scheme Hence, better results might be expected as the order is increased for problems that combine discontinuities and smooth profiles. In Fig. 3.9, similar results are shown when using the time integration scheme proposed in Sect. 3.4. A small degradation of the results can be seen in Fig. 3.9(a) as we increase the discretization order. In a similar trend, we observe less improvement in Fig. 3.9(b) than in Fig. 3.8(b). We attribute this degradation to the time partitions, which becomes more evident as the subdomains are smaller. In particular, at the boundary of each partition the method might slightly increase the amount of diffusion introduced. At these boundaries, the shock detector rely on a smaller domain to determine if the DMP is satisfied. Therefore, it is more likely to introduce more diffusion. On the other hand, the partition itself modifies the scheme introducing some error as shown in Sect. 3.6.1. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u 1st 2nd 3rd 4th IC (a) Solutions increasing the order while keeping fixed the number of control points. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u 1st 2nd 3rd 4th IC (b) Solution increasing the order using the krefinement process described above. Figure 3.8: Solution of problem (3.13) at t= 0.5for first to fourth order discretizations. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u 1st 2nd 3rd 4th IC (a) Solution using fixed number of control points, and time integration defined in Sect. 3.4 with 5 partitions. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u 1st 2nd 3rd 4th IC (b) Solution using k-refinement, and time integration defined in Sect. 3.4 with 5 partitions. Figure 3.9: Solution of problem (3.13) at t= 0.5for first to fourth order discretizations.
3.6. Numerical experiments 59 3.6.5 Boundary layer In this section the effect of the discretization order in a convection–diffusion problem is analyzed. To this end, we solve a problem with the propagation of a sharp layer and a boundary layer. In particular, we solve for Ω. = [0,1]2 (−10−4∆u+v·∇u= 0 in Ω u=uDat ∂Ω,(3.14) where v= (cos θ, sin θ),θ=−π/3, and the boundary conditions are defined as uD=(1 2+1 πarctan 10−4(y−5/6)if y= 0 0otherwise . For this test, we use the following settings: a nonlinear tolerance of 10−8,q= 2, ε= 10−8,σ= 10−6, and γ= 10−10. In Fig. 3.10(a), the solution for p= 4 is depicted. The converged solution does not exhibit any oscillation. Very sharp layers are obtained for this parameter setting. In Fig. 3.10(b), we show the profile of the solution at y= 0.1 for different orders. In this case, we start with a discretization of 50 control points per direction. Then, we increase the order using the k-refinement used previously. As previously observed for transient problems, we observe an improvement of the solution as the order is increased. (a) 3D representation of the solution. 0.0 0.2 0.4 0.6 0.8 1.0 x 0.0 0.2 0.4 0.6 0.8 1.0 u 1st order 2nd order 3rd order 4th order (b) Profiles at y= 0.1. Figure 3.10: Solution of problem (3.14) using scheme (3.6), and different discretization orders. 3.6.6 Three Body rotation Finally, we solve the transient pure convection problem (3.13) in Ω×(0,1] for Ω = [0,1]2, with v= (−2π(y−0.5),2π(x−0.5)). Initial conditions are given in [56]. Its interpolation in a first order 200 ×200 control point mesh is depicted in Fig. 3.11(a). The analytical solution of this problem is simply the translation of the profiles in the direction of the
60 Chapter 3. Arbitrary order space–time monotonicity preserving scheme convection. In particular, for t= 1, one revolution is completed and the solution is equal to the initial conditions. The purpose of this test is to evaluate how diffusive is the proposed scheme. We perform this evaluation evolving the solution until t= 1 and comparing the results with the initial conditions. The solution is computed using scheme (3.6) in combination with the shock detector in (3.10). We use the following parameters for the stabilization: q= 10,σ= 10−6, ε= 10−8, and γ= 10−10. Different meshes, time partitions, and discretization orders are used in this experiment. We start with a linear discretization of 100 ×100 control points in space, and 500 in time divided in 125 subdomains. Then, we increase the discretization order to p= 2 using the k-refinement. In order to compare first and second order discretizations, but using a similar number of control points we use a discretization with 200 ×200 control points in space, and 1000 in time divided in 250 subdomains. Finally, we assess the effect of the partitions in the temporal direction. We compare the previous discretization of 100 ×100 ×500 control points divided in 125 subdomains, with the same discretization divided in 250 subdomains. We do the same comparison for the second order discretization using 125 subdomains and when it is divided in 250 subdomains. (a) Initial conditions of the 3 body rotation. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 s 0.0 0.2 0.4 0.6 0.8 1.0 uh (b) Profile at p(x−0.5)2+ (y−0.5)2= 0.25. Figure 3.11: Three body rotation test initial conditions. Fig. 3.13 shows the solutions for 100 ×100 meshes, and 125 subdomains in time, whereas Fig. 3.14 show the ones for 250 subdomains. In both cases, a great improvement can be observed as we increase the discretization order. However, the computational cost is also increased. It is interesting to compare the solutions for first and second order discretizations using meshes with similar amount of control points, namely solutions at Fig. 3.12 and 3.13(b). For this particular problem, using a higher order discretization with similar number of control points does not improve the solution, which it is actually slightly more diffusive for p= 2. It is also worth mentioning that increasing the discretization order does not modify the behavior of the solution in terms of clipping or
3.6. Numerical experiments 61 terracing. Comparing Figs. 3.14(a) and 3.13(a), we observe that the scheme becomes more dissipative as the number of partitions is increased. This is even clearer in Fig. 3.15, where the profile of the solution at s. ={(x, y) : p(x−0.5)2+ (y−0.5)2= 0.25} is depicted. (a) 3D view. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 s 0.0 0.2 0.4 0.6 0.8 1.0 uh (b) Profile at p(x−0.5)2+ (y−0.5)2= 0.25. Figure 3.12: Three body rotation test results at t= 1 using scheme (3.6), q= 10,σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 200×200×1000 control points is used with 250 subdomains in the temporal direction. (a) Solution for a first order discretization. (b) Solution after one k-refinement. Figure 3.13: Three body rotation test results at t= 1 using scheme (3.6), q= 10,σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 125 subdomains in the temporal direction have been used.
62 Chapter 3. Arbitrary order space–time monotonicity preserving scheme (a) Solution for first order discretization. (b) Solution after one k-refinement. Figure 3.14: Three body rotation test results at t= 1 using scheme (3.6), q= 10,σ= 10−6,ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 250 subdomains in the temporal direction have been used. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 s 0.0 0.2 0.4 0.6 0.8 1.0 uh 250 subdomains 125 subdomains (a) Solution for first order discretization. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 s 0.0 0.2 0.4 0.6 0.8 1.0 uh 250 subdomains 125 subdomains (b) Solution after one k-refinement. Figure 3.15: Three body rotation test profiles for t= 1 at p(x−0.5)2+ (y−0.5)2= 0.25 using scheme (3.6), q= 10,σ= 10−6, ε= 10−8, and γ= 10−10. A first order discretization of 100 ×100 ×500 control points is used. The second order discretization is obtained using k-refinement. 125 and 250 subdomains in the temporal direction have been used. 3.7 Conclusions In this chapter an extension of the stabilization in Chapter 2 to isogeometric analysis methods have been developed. The proposed method is unconditionally DMP preserving for arbitrary high-order discretizations in space and time without any CFL-like condition. Furthermore, it is shown to be linearity-preserving in a space–time sense. Moreover, the regularized version is shown to yield better convergence behavior, especially when for the hybrid Picard–Newton method.
4.2. Preliminaries 69 subject to appropriate initial conditions u(x, 0) = u0(x). Note that the double contraction is applied as f0(u) : ∇v=Pk,γ f0 k(u)βγ vγ,β. In combination with the FE spaces described above for the spatial discretization the method of lines is being applied. The solution is approximated using u≈uh= Pi∈˜ Nh,1≤β≤mϕβ iuβ i=Pi∈˜ Nhϕiui. In a similar manner, the fluxes are approximated as f≈fh=Pi∈˜ Nh,1≤β≤mϕβ if(ui)β=Pi∈˜ Nhϕif(ui). For the sake of brevity, we use Backward Euler (BE) for the temporal discretization. Higher order time discretizations can be achieved using SSP RK methods (see [37]). In the latter case, a CFL-like condition arises to ensure that monotonicity properties in Sect. 4.2.3 are satisfied [60, 67]. The semi-discrete Galerkin FE approximation of problem (4.2) reads: find uh∈Vh such that uβ h= ¯uβ hon Γβ in,uh=u0hat t= 0, and (∂tuh,vh)+(uh,f0 h(uh) : ∇vh)−(uh,nΓout ·f0 h(uh)vh)Γout = (g,vh),∀vh∈Vh0, where uβ hand u0hare admissible FE approximations of uβand u0. In this context, we consider admissible any approximation that satisfies the maximum principle, i.e. it does not introduce new extrema. To obtain the fully discrete problem, we consider a partition of the time domain (0, T]into nts sub-intervals of length (tn, tn+1]. Then, at every time step n= 0, . . . , nts −1, the discrete problem consists in solving MδtUn+1 +KUn+1 =G, where Un+1 . = [un+1 1, ..., un+1 N]Tis the vector of nodal values at time tn+1,δt(U). = ∆t−1 n+1(Un+1 −Un), and ∆tn+1 . = (tn+1 −tn). The m×m-matrices relating nodes i, j ∈ Nhare given by Mβγ ij . = (ϕj, ϕi)δβγ, Kβγ ij . = (ϕjδβξ,f0 k(un+1 j)ξη ·∂kϕiδηγ)−(ϕjδβξ, nk·f0 k(un+1 j)ξηϕiδηγ)Γout , Gβ i . = (gβ, ϕi), where Einstein summation applies, β, γ, ξ, η ∈ {1, . . . , m}are the component indices, and δβγ is the Kronecker delta. 4.2.3 Stabilization properties In this section, we introduce some concepts required for discussing the stabilization method presented in subsequent sections. In the case of hyperbolic systems of equations, some stabilization methods are based on schemes developed for scalar equations. Let us recall some definitions used for scalar problems. Definition 4.2.1 (Local Discrete Extremum).The function vh∈Vhhas a local discrete minimum (resp. maximum) on i∈ Nhif ui≤uj(resp. ui≥uj)∀j∈ Nh(Ωi).
70 Chapter 4. Local bounds preserving FEs for first order conservation laws Definition 4.2.2 (Local DMP).A solution uh∈Vhsatisfies the local discrete maximum principle if for every i∈ Nh min j∈Nh(Ωi)\{i}uj≤ui≤max j∈Nh(Ωi)\{i}uj. Definition 4.2.3 (LED).A scheme is local extremum diminishing if, for every uithat is a local discrete maximum (resp. minimum), dui dt≤0,resp. dui dt≥0, is satisfied. One possible strategy to satisfy the above properties consist on designing a scheme that yields a positive diagonal mass matrix and a stiffness matrix that satisfies X j Aij = 0,and Aij ≤0i6=j. (4.3) In this case, it is possible to rewrite the system as mi dui dt +X j∈Nh(Ωi)\{i} Aij(uj−ui)=0,∀i∈ Nh.(4.4) As shown in [28] and [67], such a scheme satisfies the local DMP for steady problems and it is also LED when applied to transient problems. The extension of these properties to hyperbolic systems is based on analyzing them in characteristic variables. Let us consider a one-dimensional linear hyperbolic system with a constant Jacobian flux, f0. In this particular case, the continuous system can be diagonalized. Thus it is possible to discretize and solve for the characteristic variables. For example, for the set of characteristics variables, say W, the continuous system reads: ∂tW+ Λ∂xW= 0,(4.5) where Λ = diag(λ1, ..., λm)is a diagonal mby mmatrix. At this point, one can see the system as a set of independent scalar transport problems. Thus, it leads to a system with diagonal blocks after discretizing it with FEs. Assuming conditions (4.3) are satisfied for every component of problem (4.5), then the scheme will be LED for each characteristic variable. Notice that this is equivalent to forcing the original (coupled) FE approximation to have negative semi-definite offdiagonal blocks. That is, the FE discretization of the problem in characteristic variables reads (ϕj, ϕi)∂tWj+ Λ(∂xϕj, ϕi)Wj= 0.(4.6) Since in this case it is a one dimensional linear problem, we can recover the original
4.2. Preliminaries 71 problem using the fact that W=R−1U, and f0=RΛR−1. Multiplying (4.6) at the left by R, (ϕj, ϕi)∂tRR−1Uj+RΛ(∂xϕj, ϕi)R−1Uj= 0. In this case, (∂xϕj, ϕi)is simply a scalar value. Hence, we are able to recover the original (coupled) problem FE discretization. (ϕj, ϕi)∂tUj+f0(∂xϕj, ϕi)Uj= 0. Thus, if f0(∂xϕj, ϕi)is negative semi-definite for j6=i, then the problem in characteristic variables will satisfy conditions (4.3) for each variable. In the case of more general multidimensional problems (e.g. Euler equations), this would only imply that the scheme is LED for a certain set of local characteristic variables. Furthermore, if the flux Jacobian f0is not linear, then even the definition of the matrix Aij (relating nodes iand j) is not trivial. Let us recall the definition of these blocks for Euler equations Kij . = (ϕj,f0 k(uj)·∂kϕi)−(ϕj, nk·f0 k(uj)ϕi)Γout =−f0 k(uj)·(∂kϕj, ϕi), where we have undone integration by parts. It is easy to check that Pj(∂kϕj, ϕi)=0. Hence, we can write X j∈Nh(Ωi) Kijuj=X j∈Nh(Ωi)\{i}−(∂kϕj, ϕi)(f0 k(uj)·uj−f0 k(ui)·ui). As previously stated, it is not straightforward in the case of Euler equations to rewrite the discrete problem in the form of (4.4). However, making use of special density-averaged variables it is possible to rewrite the previous expression as X j∈Nh(Ωi) Kijuj=X j∈Nh(Ωi)\{i}−f0 k(uij)·(∂kϕj, ϕi)(uj−ui), where uij are the Roe mean values [82]. For an ideal gas, these are defined as ρij =√ρiρj,mij =mi√ρj+mj√ρi √ρi+√ρj ,(ρE)ij =1 2−γρijHij −|mij|2 2ρij , where Hij is the average enthalpy Hij =Hi√ρi+Hj√ρj √ρi+√ρj ,and Hi=E−pi ρi . Therefore, using this density-averaged variables it is possible to rewrite Euler problem in the form of (4.4). Hence, if −f0 k(uij)·(∂kϕj, ϕi)has non-positive eigenvalues, then the scheme will be LED for a certain set of local characteristic variables. Schemes that satisfy this property are named local bounds preserving schemes in the literature
72 Chapter 4. Local bounds preserving FEs for first order conservation laws [75]. This reasoning above motivated the definition of the LED principle for hyperbolic systems of equations by Kuzmin [59] and coworkers. Adapted from this principle, we define local bounds preserving schemes as follows. Definition 4.2.4. The semi-discrete scheme X j Mij∂tuj+X j6=i Aij(uj−ui) = 0 is said to be local bounds preserving if Mis diagonal with positive entries (i.e. Mij = miδijIm×m), Aij has non-positive eigenvalues for every j6=i, and PjAij =0. Unfortunately, to the best of our knowledge, satisfying this definition does not ensure positivity of density, internal energy, or non-decreasing entropy. In any case, numerical schemes based on this definition have shown good numerical behavior [59, 62, 70, 75]. Several stabilization strategies have been defined based on the previous ideas. One of the most simple strategies consists on adding a scalar artificial diffusion term proportional to the spectral radius of Aij [59, 73]. Sometimes, this strategy is named Rusanov artificial diffusion, since for linear FEs in one dimension the scheme results in the Rusanov Riemann solver [59, 91]. Without any special treatment, the resulting scheme is only first order accurate. The key to recovering high-order convergence is to modulate the action of the artificial diffusion term, and restrict its action to the vicinity of discontinuities. We base our stabilization term in Rusanov artificial diffusion and a differentiable shock detector recently developed for scalar problems in Chapters 2 and 3. 4.3 Nonlinear stabilization As previously discussed, the Galerkin FE discretization yields oscillatory solutions in regions around discontinuities. We supplement the original scheme with an artificial diffusion term to stabilize it and mitigate these oscillations. The proposed stabilization term is given by Bh(wh;uh,vh). =X Ke∈ThX i,j∈Nh(Ke) 1≤β,γ≤m νe ij(wh)`(i, j)vβ i·δβγuγ j,(4.7) for any uh∈Vhand vh∈Vh0. Here, `(i, j). = 2δij −1is a graph Laplacian operator defined in Chapter 2, and νe ij(wh)is the element-wise artificial diffusion defined as νe ij(wh). = max αi(wh)λmax ij ,αj(wh)λmax ji ,for j∈ Nh(Ωi)\{i}, νe ii(wh). =X j∈Nh(Ωi)\{i} νe ij(wh),(4.8) where λmax ij is the spectral radius of the elemental convection matrix relating nodes i, j ∈ Nh, i.e. ρf0(uij)·(∇ϕj, ϕi)Ke. As previously introduced, this artificial diffusion
4.3. Nonlinear stabilization 73 term is based on Rusanov scalar diffusion [61]. It is important to mention that the eigenvalues of these matrices can be easily computed as λ1,..,d =vij ·ce ij, λd+1 =vij ·ce ij −ckce ijk, λd+2 =vij ·ce ij +ckce ijk(4.9) where ce ij = (∇ϕj, ϕi)Ke,and c=v u u t(γ−1) Hij −kmijk2 2ρ2 ij !. We denote by αi(wh)the shock detector used for modulating the action of the artificial diffusion term. The idea behind the definition of this detector is minimizing the amount of artificial diffusion introduced while stabilizing any oscillatory behavior. In regions where the local DMP (see Def. 4.2.2) is not satisfied for any chosen set of components, we ensure that Def. 4.2.4 is satisfied. αi(wh)must be a positive real number which takes value 1 when uh(xi)is an inadmissible value of uh, and smaller than 1 otherwise. To this end, we define αi(uh). = max{αi(uβ h)}β∈C,(4.10) where Cis the set of components that are used to detect inadmissible values of uh, e.g. density and total energy in the case of Euler equations. For simplicity, we restrict ourselves to the components of uh. However, derived quantities such as the pressure or internal energy can be also used. In order to introduce the shock detector, let us recall some useful notation from Chapter 2. Let rij =xj−xibe the vector pointing from node xito xjwith i, j ∈ Nhand ˆ rij . =rij |rij |. Recall that the set of points xjfor j∈ Nh(Ωi)\{i}define the macroelement Ωiaround node xi. Let xsym ij be the point at the intersection between ∂Ωiand the line that passes through xiand xjthat is not xj(see Fig. 4.1). The set of all xsym ij for all j∈ Nh(Ωi)\{i}is represented with Nsym h(Ωi). We define rsym ij . =xsym ij −xi. Given xsym ij in two dimensions, let as call aand bthe indices of the vertices such that they define the edge in ∂Ωithat contains xsym ij . We define usym jas the value of uhat xsym ij , i.e. uh(xsym ij ). Figure 4.1: usym drawing. Both usym ij and xsym ij are only required to construct a linearity preserving shock detector. Let us define the jump and the mean of a linear approximation of component β
74 Chapter 4. Local bounds preserving FEs for first order conservation laws of the unknown gradient at node xiin direction rij as r∇uβ hzij . =uβ j−uβ i |rij|+usym,β j−uβ i |rsym ij |, nn|∇uβ h·ˆ rij|ooij . =1 2 |uβ j−uβ i| |rij|+|usym,β j−uβ i| |rsym ij |!. For each component in C, we use the same shock detector developed in Chapter 2. Let us recall its definition αi(uβ h). = Pj∈Nh(Ωi)r∇uβ hzij Pj∈Nh(Ωi)2nn∇uβ h·ˆ rijooij q if Pj∈Nh(Ωi){{|∇·ˆ rij|}}ij 6= 0 0otherwise . (4.11) From Lm. 2.3.1 we know that (4.11) is valued between 0 and 1, and it is only equal to one if uβ h(xi)is a local discrete extremum (in a space–time sense as in Def. 4.2.1). Since the linear approximations of the unknown gradients are exact for uβ h∈ P1, the shock detector vanishes when the solution is linear. Thus, it is also linearly preserving for every component in C. This result follows directly from Th. 2.4.5. The final stabilized problem in matrix form reads as follows. Find uh∈Vhsuch that uh=uhon ∂Ω,uh=u0hat t= 0, and M(un+1 h)δtUn+1 +Kij(un+1 h)Un+1 =G(4.12) for n= 1, ..., nts, where Mij(un+1 h). = [1 −max (αi,αj)] (ϕj, ϕi)Im×m+ max (αi,αj) (δij, ϕi)Im×m, Kij(un+1 h). =Kij +Bij, and Bij(uh) = Bh(uh;ϕj, ϕi), for i, j ∈ Nh. Lemma 4.3.1 (Local bounds preservation).Consider uh∈Vhwith component βin the set of tracked variables C. The stabilized problem (4.12) is local bounds preserving as defined in Def. 4.2.4 at any region where uβ hhas extreme values. Proof. If component β∈Cof uhhas an extremum at xi, then from Lm. 2.3.1 we know that αi(uβ h) = 1. Moreover, from (4.10) is easy to see that αi(uh) = 1. In this case, Mij(uh)=(δij, ϕi)Im×m. Hence, Mij(uh)=0for j6=iand Mii(uh) = mi. Therefore, we can rewrite the system as follows mi∂tui+X j∈Nh(Ωi)\{i} Kij(uij)(uj−ui) = mi∂tui+X j∈Nh(Ωi)\{i}X Ke∈Thf0(uij)·(∇ϕj, ϕi)Ke−νe ijIm×m(uj−ui) = 0.
4.3. Nonlinear stabilization 75 We need to prove that the eigenvalues of Kij(uij)are non-positive. To this end, let us show the following inequality holds X Ke∈Th ρf0(uij)·(∇ϕj, ϕi)Ke≥ρ(f0(uij)·(∇ϕj, ϕi)). From (4.9), it is easy to check that ρf0(uij)·(∇ϕj, ϕi)Ke=|vij ·ce ij|+cijkce ijk. Since cij =PKe∈Thce ij, we have that X Ke∈Thvij ·ce ij≥ |vij ·cij|,and X Ke∈Th cijkce ijk ≥ cijkcijk. Hence, Peρ(Ke ij(uij)) ≥ρ(Kij(uij)). Moreover, by definition (see (4.8)), νe ij ≥ρf0(uij)·(∇ϕj, ϕi)Kefor j6=i. Furthermore, from (4.7), is easy to see that ρ(Be ij(uij)) ≥ρ(Ke ij(uij)). Therefore, ρ(Bij(uij)) ≥ρ(Kij(uij)). Finally, since Kij =Kij +Bij and Bij =PeBe ij = Pe−νe ijIm×mfor all j6=i, then the maximum eigenvalue of Kij(uij)is non-positive, which completes the proof. 4.3.1 Differentiability In the case of steady, or implicit time integration, differentiability plays a role in the convergence behavior of the nonlinear solver. This is especially important if one wants to use Newton’s method. In the case of scalar problems it has been shown in the previous chapters and [5] that convergence is greatly improved after few modifications to make a scheme twice-differentiable. In this section, we introduce a set of regularizations applied to all non-differentiable functions present in the stabilized scheme introduced above. In order to regularize these functions, we follow a similar strategy as in Chapter 2. Absolute values are substituted by |x|1,εh=px2+εh,|x|2,εh=x2 px2+εh . Note that |x|2,εh≤ |x|≤|x|1,εh. Next, we also use a smooth maximum function, max σh(·), as max σh(x, y). =|x−y|1,σh 2+x+y 2≥max(x, y).(4.13) In addition, we need a smooth function to limit the value of any given quantity to one. To this end, we use Z(x). =(2x4−5x3+ 3x2+x, x < 1, 1, x ≥1.
76 Chapter 4. Local bounds preserving FEs for first order conservation laws The set of twice-differentiable functions defined above allow us to redefine the stabilization term introduced in Sect. 4.3. In particular, we define ˜ Bh(wh;uh,vh). =X Ke∈ThX i,j∈Nh(Ke) 1≤β,γ≤m ˜νe ij(wh)`(i, j)vβ i·δβγuγ j, where ˜νe ij(wh). = max σhαεh,i(wh)λmax ij ,αεh,j(wh)λmax ji ,for j∈ Nh(Ωi)\{i}, ˜νe ii(wh). =X j∈Nh(Ωi)\{i} ˜νe ij(wh).(4.14) Let us note that λmax ij needs to be regularized as λmax ij =vijce ij1,εh +ckce ijk. The shock detector is also redefined to use the regularized version of the shock detector, which reads αεh,i(uh). = max σh{αεh,i(uβ h)}β∈C. In the case of the component shock detector we recall the definition in Chapter 2 αεh,i(uβ h). = Z Pj∈Nh(Ωi)r∇uβ hzij1,εh +ζh Pj∈Nh(Ωi)2∇uβ h·ˆ rij2,εhij +ζh q ,(4.15) where ζhis a small value for preventing division by zero. Finally, the twice-differentiable stabilized scheme reads: Find uh∈Vhsuch that uh=uhon ∂Ω,uh=u0hat t= 0, and ˜ M(un+1 h)δtUn+1 +˜ Kij(un+1 h)Un+1 =Gfor n= 1, ..., nts,(4.16) where ˜ Mij(un+1 h). = [1 −max σh(αεh,i,αεh,j)] (ϕj, ϕi)Im×m + max σh(αεh,i,αεh,j) (δij, ϕi)Im×m, ˜ Kij(un+1 h). =Kij(un+1 h) + ˜ Bij(un+1 h), and ˜ Bij(uh) = ˜ Bh(uh;ϕj, ϕi), for i, j ∈ Nh. Corollary 4.3.2. The differentiable scheme in Eq. (4.14) is local bounds preserving, as defined in Def. 4.2.4, at any region where uβ hhas extreme values for every βin C. Proof. For an extreme value of uβ h, since |x|2,εh≤ |x|≤|x|1,εhthe quotient of (4.15) is larger than one. Hence, by definition of Z(x),αεh,i is equal to 1. At this point, it is easy to check that ˜νe ij ≥νe ij in virtue of the definition of max σh. Therefore, ρ(˜ Be ij(uh)) ≥ ρ(Be ij(uh)), completing the proof.
4.4. Nonlinear solver 77 Moreover, it is important to mention that the differentiable shock detector is weakly linearly-preserving as ζhtends to zero. This result follows directly from Sect. 2.7. In order to obtain a differentiable operator, we have added a set of regularizations that rely on different parameters, e.g., σh, εh, ζh. Giving a proper scaling of these parameters is essential to recover theoretic convergence rates. In particular, we use the following relations σh=σ|λmax|2L2(d−3)h4, εh=εL−4h2, ζh=L−1ζ, (4.17) where dis the spatial dimension of the problem, Lis a characteristic length, and σ, ε, and ζare of the order of the unknown. 4.4 Nonlinear solver In this section, we describe the method used for solving the nonlinear system of equations arising from the scheme introduced above. In particular, we use a hybrid Picard–Newton approach in order to increase the robustness of the nonlinear solver. Moreover, for the differentiable version we also use a continuation method to improve the nonlinear convergence. We represent the residual of the equation (4.16) at the k-th iteration by R(uk,n+1 h), i.e., R(uk,n+1 h). =˜ M(uk,n+1 h)δtUk,n+1 +˜ Kij(uk,n+1 h)Uk,n+1 −G.(4.18) Hence, the Jacobian is defined as J(uk,n+1 h). =∂R(uk,n+1 h) ∂Uk,n+1 (4.19) =˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h) + ∂˜ M(uk,n+1 h) ∂Uk,n+1 δtUk,n+1 +∂˜ Kij(uk,n+1 h) ∂Uk,n+1 Uk,n+1. Therefore, Newton method consist on solving J(uk,n+1 h)∆Uk+1,n+1 =−R(uk,n+1 h). However, it is well known that Newton method can diverge if the initial guess of the solution u0,n+1 his not close enough to the solution. In order to improve the robustness, we introduce the following modifications. We use a line–search method to update the solution at every time step. Thus, the new approximation is computed as Uk+1,n+1 =Uk,n+1 +λ∆Uk+1,n+1, where λis computed (approximately) such that it minimizes kR(uk+1,n+1 h)k. As introduced at the beginning of the section, we also use a hybrid approach combining Newton method with Picard linearization. Picard nonlinear iterator can be obtained removing the last two terms of (4.19), i.e., ˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h)∆Uk+1,n+1 =−R(uk,n+1 h).(4.20)
78 Chapter 4. Local bounds preserving FEs for first order conservation laws Clearly, it is equivalent to ˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h)Uk+1,n+1 =˜ M(uk,n+1 h)Un+G. Moreover, we modify the definition of left hand side terms in (4.20) to enhance the robustness of the method. In particular, we use αi= 1 for computing these terms while we use the value obtained from (4.10) for the residual. Using this strategy the solution remains unaltered, but the obtained approximations uk,n+1 hfor intermediate values of k are more diffusive. Even though this modification slows the nonlinear convergence, it is essential at the initial iterations. Otherwise, the robustness of the method might be jeopardized. The resulting iterative nonlinear solver consists in the following. We iterate using Picard method in (4.20), with the modification described above, until the L2norm of the residual is smaller than a given tolerance. In this chapter, we use tolerances close to 10−2. Afterwards, Newton method with the exact Jacobian in (4.19) is used until the desired nonlinear convergence criteria is satisfied. For the differentiable stabilization, we also equip the above scheme with a continuation method on the regularization parameters. In order to accelerate the convergence of the method, we use high values for the parameters during the first iterations. This results in a more diffusive solution, but nonlinear convergence is accelerated. As the nonlinear approximation is closer to the solution, we diminish the value of the parameters to avoid introducing excessive artificial diffusion to the system. This process is preformed gradually as a function of the residual in (4.18). In particular, we use the following relation εk= ˜εkR(uk,n+1 h)k kR(u0,n+1 h)k, where εkis the effective parameter used in relations 4.17, and ˜εis parameter defined by the user. We summarize the nonlinear solver introduced above in Alg. 3. 4.5 Numerical experiments In this section, we perform several numerical experiments to assess the numerical scheme introduced in the previous sections. First, we perform a convergence analysis to assess its implementation. Then, we use a steady benchmark test to analyze the effectiveness of the regularization parameters. We also analyze their effectiveness in the case of a transient problem. Finally, we solve a slightly more challenging steady benchmark test. In all experiments below we assume that the ideal gas state equation applies, and we use an adiabatic index of γ= 1.4. From previous experience [4, 5, 16, 18], the effects of parameters σand εto the nonlinear convergence and numerical error are analogous. Hence, we consider ε= 10−2σ. In addition, for all the tests below, the density is
4.5. Numerical experiments 85 qin Fig. 4.10(c). Therefore, one can come to the conclusion that in order to achieve a given accuracy it is preferable to use the differentiable scheme with a slightly larger value of qrather than the non-differentiable scheme and a low value for q. 4.5.4 Scramjet Finally, we solve a problem with a supersonic flow that develops a complex shock pattern. This test consists of a M= 3 channel that narrows along the streamline and has two internal obstacles. In particular, Fig. 4.11 is an illustration of the domain and Tab. 4.3 lists the coordinates of the points defining the domain. The problem is solved directly to steady state, and two different meshes have been used. The coarsest mesh used has 18476 Q1elements and the finest mesh has 63695 Q1elements. Figure 4.11: Scramjet test scheme. Table 4.3: Domain coordinates for the scramjet test. Wall a b c d e f xi0.0 0.4 4.9 12.6 14.25 16.9 yi3.5 3.5 2.9 2.12 1.92 1.7 Interior obstacle A B C D E xi4.9 8.9 9.4 12.6 14.25 yi-1.4 -0.5 -0.5 -1.4 -1.2 In order to solve this problem, the hybrid nonlinear solver described in Sect. 4.4 is used with the help of the continuation scheme. The tolerance for switching from the Picard to Newton linearization is set to 5·10−2. The nonlinear convergence criteria for this benchmark is k∆(uk+1 h)k kuk hk<10−6. We also set a maximum number of iterations of 500. In this test, we use q={2,5},γ= 10−10,˜ε={1,10−2,10−4}, and εk=σk10−2. Even though σ= 102might seem a high value, we recall that it is used in the context of a continuation method. Therefore, the effective value of σkis lower than 1 for the converged solution. Moreover, the actual value used in (4.13) is computed using the relations in (4.17). Figs. 4.12–4.13 show, respectively, the Mach and density contours for the fine mesh, q= 5,˜ε= 1,˜ε= 102, and γ= 10−10. The nonlinear convergence history for this
86 Chapter 4. Local bounds preserving FEs for first order conservation laws configuration is depicted in Figs. 4.17(c)–4.17(d). The obtained values for the Mach number and the density are comparable to those in [61, 75]. The shocks are well resolved. Even when using q= 2 the shocks are properly resolved and only slightly more smeared than for q= 5, see Fig. 4.14. If instead, the coarse mesh is used (see Fig. 4.15), the solution is more dissipative. However, the scheme is able to capture most of the features present in the solution. Figure 4.12: Scramjet Mach contours when a mesh of 63695 Q1ele- ments is used, with parameters q= 5,γ= 10−10, and ˜ε= 1. Figure 4.13: Scramjet Mach contours when a mesh of 63695 Q1ele- ments is used, with parameters q= 5,γ= 10−10, and ˜ε= 1. Figure 4.15: Scramjet Mach contours when a mesh of 18476 Q1ele- ments is used, with parameters q= 2,γ= 10−10, and ˜ε= 1.
4.5. Numerical experiments 87 Figure 4.14: Scramjet Mach contours when a mesh of 63695 Q1ele- ments is used, with parameters q= 2,γ= 10−10, and ˜ε= 1. Figs. 4.16–4.17 show the nonlinear convergence history in terms of the relative residual reduction and the relative solution increment between iterations. We can observe that the convergence is not ensured for an arbitrary choice of the regularization parameters. In fact, only the tests that use ˜ε= 1 do not diverge for q= 5 , regardless of the mesh used. Therefore, we can see that increasing the values of the regularization parameters not only improves the convergence, but also the robustness of the method. 0 100 200 300 400 500 iteration 10−3 10−2 10−1 100 kR(uk)k/kR(u0)k ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (a) q= 2. 0 100 200 300 400 500 iteration 10−5 10−4 10−3 10−2 10−1 100 kuk+1 −ukk/kukk ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (b) q= 2. 0 100 200 300 400 500 iteration 10−3 10−2 10−1 100 101 kR(uk)k/kR(u0)k ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (c) q= 5. 0 100 200 300 400 500 iteration 10−5 10−4 10−3 10−2 10−1 100 kuk+1 −ukk/kukk ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (d) q= 5. Figure 4.16: Comparison of the convergence behavior for the Scramjet test and different regularization parameters choices. A coarse mesh of 18476 Q1elements is used.
88 Chapter 4. Local bounds preserving FEs for first order conservation laws However, it is important to mention that even if we can improve the convergence behavior of these types of methods, this is not enough for directly solving to steady state problems with complex shock patterns. For instance, even if the solution Fig. 4.15 seems to be correct, the scheme was unable to converge to the desired tolerance (see Figs. 4.16(c) and 4.16(d)). However the ability to introduce differentiability into the definition of the shock detector, for robustness and increased nonlinear convergence rates, could be coupled with popular pseudo-time stepping approaches [51, 86] to pursue improved methods for complex shock type systems. 0 100 200 300 400 500 iteration 10−3 10−2 10−1 100 kR(uk)k/kR(u0)k ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (a) q= 2. 0 100 200 300 400 500 iteration 10−6 10−5 10−4 10−3 10−2 10−1 100 kuk+1 −ukk/kukk ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (b) q= 2. 0 50 100 150 200 250 300 350 400 450 iteration 10−4 10−3 10−2 10−1 100 kR(uk)k/kR(u0)k ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (c) q= 5. 0 50 100 150 200 250 300 350 400 450 iteration 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 kuk+1 −ukk/kukk ˜ε= 0 ˜ε= 10−4 ˜ε= 10−2 ˜ε= 1 (d) q= 5. Figure 4.17: Comparison of the convergence behavior for the Scramjet test and different regularization parameters choices. A fine mesh of 63695 Q1elements is used. 4.6 Conclusions In this chapter, a differentiable local bounds preserving stabilization method for Euler equations has been developed. This stabilization is based on the combination of a differentiable shock detector, a partially lumped mass matrix, and Rusanov artificial diffusion operator. The resulting scheme has been successfully tested for steady and transient benchmark problems. Numerical results show that the proposed method exhibits good stability
4.6. Conclusions 89 properties. Furthermore, it is able to provide well resolved sharp shocks in both steady and transient problems. In addition, to improve nonlinear convergence, a continuation method for the regularization parameters present in the differentiable stabilization has also been proposed. Nonlinear convergence of the scheme has been analyzed for the differentiable version and compared with its non-regularized counterpart. In general terms, the differentiable stabilization shows better convergence, especially when the hybrid Picard–Newton method is used. For small steady problems, the scheme is able to converge directly to the steady state solution without making use of pseudo-transient time stepping. However, for problems with complex shock patterns the scheme only converges to moderate tolerances. Numerical results also show that differentiability not only can improve nonlinear convergence, but it also improves the robustness of the method. In the case of transient problems, some improvement in the computational cost is observed. However, since the non-differentiable method already exhibits good nonlinear convergence, there is not much room for improvement. Nevertheless, it is possible to show that the differentiable stabilization can achieve a similar accuracy while requiring a lower computational cost.
Chapter 5 Monotonicity-preserving FE schemes with AMR for hyperbolic problems This chapter is focused on the extension and assessment of the monotonicity-preserving scheme in Chapter 2 and the local bounds preserving scheme in Chapter 4 to hierarchical octree AMR. Whereas the former can readily be used on this kind of meshes, the latter requires some modifications. A key question that we want to answer in this chapter is whether to move from a linear to a nonlinear stabilization mechanism pays the price when combined with shock-adapted meshes. Whereas nonlinear (or shockcapturing) stabilization leads to improved accuracy compared to linear schemes, it also negatively hinders nonlinear convergence, increasing computational cost. We compare linear and nonlinear schemes in terms of the required computational time versus accuracy for several steady benchmark problems. Numerical results indicate that, in general, nonlinear schemes can be cost-effective for sufficiently refined meshes. Besides, it is also observed that it is better to refine further around shocks rather than using sharper shock capturing terms, which usually yield stiffer nonlinear problems. In addition, a new refinement criterion has been proposed. The proposed criterion is based on the graph Laplacian used in the definition of the stabilization method. Numerical results show that this shock detector performs better than the well-known Kelly estimator for problems with shocks or discontinuities. 5.1 Introduction Natural phenomena can develop shock waves in different scenarios. A classical example is the shock wave generated by an object traveling faster than sound. The numerical modeling of problems with shocks is still a challenge, especially when the admissible physical solution has some physical constraints, e.g., positivity or non-negativity, that must be preserved at the discrete level to have well-posedness; E.g., the fluid density and temperature are positive quantities in a compressible flow. 91
92 Chapter 5. Monotonicity-preserving FE schemes with AMR Several numerical schemes have been proposed so far to approximate this kind of problems by combining FVM or dG FEs for space discretization with explicit time integrators (see [27, 34, 68, 91]). Explicit time integrators are only stable under a CFL restriction over the time step size, which implies to capture all time scales. Thus, explicit methods are not suitable for problems in which the smallest time scales are not of interest. For instance, the fastest time scales at a confined plasma in a nuclear fusion reactor are not of engineering interest whereas explicit time integration is unaffordable in practical simulations [54]. Implicit monotonicity-preserving (or at least positivity-preserving) methods are still scarce. As proved by Godunov [36], linear monotonicity-preserving schemes can be at most first-order accurate. For scalar problems (and under some mesh restrictions), Burman and Ern [25], Barrenechea and co-workers [13, 14], Kuzmin and co-workers [58, 61, 71], and Badia and Hierro [7, 8] have proposed nonlinear schemes that preserve monotonicity and can presumably attain higher order accuracy.1However, these properties comes at the cost of solving a very stiff nonlinear problem [57]. The authors [4, 5] have proposed differentiable schemes that improve the nonlinear convergence behavior of previous methods. For hyperbolic systems of equations, numerical methods are even less well developed. For explicit time integration, Guermond and Popov [42] have recently proposed a cG FE scheme that preserves positivity of density and energy under certain CFL-like condition. Unfortunately, these ideas cannot be easily extended to implicit time integration and we are not aware of any implicit method that theoretically satisfies such properties. Kuzmin and co-workers [61, 69, 74, 75] have proposed various schemes based on FCT [72] that are experimentally robust, but lack of a theoretical analysis. Besides, this strategy also yields very stiff nonlinear problems. Differentiable schemes for compressible flows have been proposed in Chapter 4 to alleviate (but not eliminate) this problem. Shocks are non-smooth and localized and thus suitable for AMR [31, 92]. AMR allows one to increase the mesh resolution only in the vicinity of shocks or discontinuities. In brief, the AMR process can be divided into two main ingredients. On the one hand, to estimate the error at each element. On the other hand, to decide which elements need to be refined or coarsened. This iterative process provides a mesh locally adapted to the features of the problem at hand. As a result, it is a nonlinear approximation scheme which tries to minimize the error for a target computational cost. If performed optimally, AMR exhibits exponential convergence even for solutions with limited regularity [31]. In this context, a key question is whether it is computationally more effective to consider a nonlinear high-order scheme (with the nonlinear convergence issues) or a cheaper linear (first-order) scheme in a much refined mesh. The motivation of this chapter is to shed light on this issue. First, we adapt the schemes developed in Chapters 2 and 4 to hierarchical octree AMR [10, 89]. Next, we propose a refinement criterion that 1In this chapter, schemes with nonlinear stabilization are also referred to as high-order and linear stabilization schemes as low or first-order.
5.2. Preliminaries 93 relies on information already present in the stabilization technique; nonlinear stabilization methods include a shock detector to activate the artificial diffusion only close to discontinuities. We propose to use a modification of the shock detector in Chapter 2 to drive the AMR process. This chapter is structured as follows. First, we introduce the problem, its discretization, and monotonicity properties for scalar problems and hyperbolic systems in Sect. 5.2. Then, the stabilization techniques are introduced in Sect. 5.3. Sect. 5.4 is devoted to the AMR strategy. We introduce the nonlinear solvers in Sect. 5.5. Finally, we show numerical experiments in Sect. 5.6 and draw some conclusions in Sect. 5.7. 5.2 Preliminaries 5.2.1 Continuous problem Let us consider an open bounded and connected domain, Ω∈Rd, where dis the number of spatial dimensions. Let ∂Ωbe the Lipschitz continuous boundary of Ω. The conservative form of a first order hyperbolic problem reads ∂tu−∇·f(u) = g,in Ω×(0, T], uβ(x, t) = ¯uβ(x, t),on Γβ in ×(0, T], β = 1, ..., m, u(x, 0) = u0(x), x ∈Ω, (5.1) where u={uβ}m β=1 are m≥1conserved variables, fis the physical flux, ¯uβ(x, t)are the boundary values for the βth-component of u,u0(x)are the initial conditions, and g(x, t)is a function defining the body forces. Note that the flux, f:Rm→Rm×d, is composed of f={fi}d i=1, where fi:Rm→Rmis the flux in the ith spatial direction. We denote by f0:Rm→Rm×m×dthe flux Jacobian. Let n∈Rdbe any direction vector. Since the system is hyperbolic, the flux Jacobian in any direction is diagonalizable and has only real eigenvalues, i.e., f0(u)·n=Pd i=1 f0 i(u)niis diagonalizable with real eigenvalues {λβ}m β=1. These eigenvalues might have different multiplicities and different signs. Hence, for a given direction n, each characteristic variable might be convected forward (along n) or backwards (along −n). Therefore, it is convenient to define inflow and outflow boundaries for each component. The inflow boundary for component βis defined as Γβ in . ={x∈∂Ω : λβ(f0(u)·n∂Ω)≤0}, where n∂Ωis the unit outward normal to the boundary and λβis the βth-eigenvalue of the flux Jacobian. We define the outflow boundary as Γβ out . =∂Ω\Γβ in. We refer the reader to [34, 43, 91] for a detailed discussion on boundary conditions for hyperbolic problems. In the present study, we will also consider the steady counterpart of (5.1), which is obtained by dropping the time derivative term and the initial conditions. In this chapter, we work with both scalar convection equations and Euler equations. Taking m= 1 and f(u). =vuwith va divergence-free convection field, we recover the well known scalar transport problem. On the other hand, Euler equations for ideal gases
94 Chapter 5. Monotonicity-preserving FE schemes with AMR are recovered by defining m=d+ 2 and u. = ρ m ρE ,f. = m m⊗v+pI v(ρE +p) ,and g. = 0 b b·v+r , where ρis the density, Eis the total energy, pis the pressure, m={m1, . . . , md}, where mi=ρvi, is the momentum, v={v1, . . . , vd}is the velocity, b={b1, . . . , bd}are the body forces, ris an energy source term per unit mass, and Iis the identity matrix of dimension d×d. In addition, the system is equipped with the ideal gas equation of state p= (γ−1)ρı, where ı=E−1 2kvk2is the internal energy and γis the adiabatic index. 5.2.2 Discretization The discretization used in this chapter is able to adapt its local size to the features of the problem at hand. In particular, it is a hierarchically refined octree-based hexahedral mesh [89]. This type of discretizations are constructed hierarchically. At every step of the refinement process, marked cells are refined into four (eight) cells in 2D (3D). The adaptation of the mesh to the problem at hand is achieved by only marking for refining a targeted amount of cells. This results in a mesh with different refinement levels at different regions. Hanging nodes appear at the interface between cells at different refinement levels. These are nodes that only belong to the cells at a higher refinement level (see Fig. 5.1). In our case, the meshes used are 2:1 balanced. This restriction implies that there can only be a difference of one refinement level between neighboring cells. This restriction is a trade-off between implementation complexity and performance gain that has been adopted by many AMR codes [89]. Figure 5.1: Example of a mesh with hanging nodes. Hanging nodes need to be treated carefully in the case of working with conforming FE discretizations. Otherwise, associating a regular degree of freedom (DOF) to a hanging node may lead to discontinuities in the approximated solution. To preserve continuity of the FE space, hanging DOFs values are not included in the assembled system of equations but obtained by interpolating the values of the neighboring regular DOFs. For more details in the definitions of these constraints we refer the reader to [9–11].
5.3. Nonlinear stabilization 101 Corollary 5.3.1 (DMP).The solution of the discrete problem (5.10) with m= 1 and using the shock detector (5.9) satisfies the local DMP in Def. 5.2.2 if g= 0 and, for every control point i∈ Nhsuch that uiis a local discrete extremum, it holds: Kij(uh)≤0,∀j∈ Nh(Ωi)\{i},X j∈Nh(Ωi) Kij(uh) = 0. Moreover, the resulting scheme is linearity-preserving as defined in Def. 5.2.5, i.e., Bij(uh) = 0 for uh∈ P1(Ωi). Proof. The stabilization scheme for scalar problems is defined on the assembled system. Hence, the modifications introduced in the assembly procedure do not affect the reasoning in the proof of [5, Thm. 5.2]. Lemma 5.3.2 (Local bounds preservation).Consider uh∈Vhwith component βin the set of tracked variables C. The stabilized problem (5.10) is local bounds preserving as defined in Def. 5.2.4 at any region where uβ hhas extreme values. Proof. If component β∈Cof uhhas an extremum at xi, we know from Lm. 2.3.1 that αi(uβ h)=1. Moreover, it is easy to see from (5.8) that αi(uh)=1. In this case, Mij(uh) = δij PjMij. Hence, Mij(uh)=0for j6=iand Mii(uh) = mi. Therefore, we can rewrite the system as follows mi∂tui+X j∈Nh(Ωi)\{i} Kij(uij)(uj−ui) = mi∂tui+X j∈Nh(Ωi)\{i}X Ke∈Th (f0(uij)) ·(∇ϕj, ϕi)Ke +X k∈M(j) Ckj(∇ϕk, ϕi)Ke +X k∈M(i) Cki(∇ϕj, ϕk)Ke +X k∈M(i)∩M(j) CkiCkj(∇ϕk, ϕk)Ke −νe ijIm×m(uj−ui) = 0.
102 Chapter 5. Monotonicity-preserving FE schemes with AMR We need to prove that the eigenvalues of Kij(uij)are non-positive. To this end, let us show that the following inequality holds X Ke∈Thρf0(uij)·(∇ϕj, ϕi)Ke +X k∈M(j) Ckj ρf0(uij)·(∇ϕk, ϕi)Ke +X k∈M(i) Cki ρf0(uij)·(∇ϕj, ϕk)Ke +X k∈M(i)∩M(j) CkiCkj ρf0(uij)·(∇ϕk, ϕk)Ke≥ρ(f0(uij)·(∇ϕj, ϕi)). From (5.6), it is easy to check that ρf0(uij)·(∇ϕj, ϕi)Ke=|vij ·ce ij|+cijkce ijk. We have that cij = (∇ϕj, ϕi) = X Ke∈Thce ij +X k∈M(j) Ckjce ik +X k∈M(i) Ckice kj +X k∈M(i)∩M(j) CkiCkjce kk, where ce ij = (∇ϕj, ϕi)Ke. Thus, X Ke∈Thvij ·ce ij+X k∈M(j) Ckj |vij ·ce ik| +X k∈M(i) Cki vij ·ce kj +X k∈M(i)∩M(j) CkiCkj |vij ·ce kk|≥ |vij ·cij|, and X Ke∈Thcijkce ijk+X k∈M(j) Ckj |vij ·cij|kce ikk +X k∈M(i) Cki |vij ·cij|kce kjk +X k∈M(i)∩M(j) CkiCkj |vij ·cij|kce kkk≥cijkcijk.
5.3. Nonlinear stabilization 103 Therefore, Peρ(Ke ij(uij)) ≥ρ(Kij(uij)). Moreover, by definition (see (5.5)), νe ij ≥ρf0(uij)·(∇ϕj, ϕi)Ke +X k∈M(j) Ckj ρf0(uij)·(∇ϕk, ϕi)Ke +X k∈M(i) Cki ρf0(uij)·(∇ϕj, ϕk)Ke +X k∈M(i)∩M(j) CkiCkj ρf0(uij)·(∇ϕk, ϕk)Kefor j6=i. Furthermore, it is easy to infer from (5.4) that ρ(Be ij(uij)) ≥ρ(Ke ij(uij)). Hence, ρ(Bij(uij)) ≥ρ(Kij(uij)). Finally, since Kij =Kij+Bij and Bij =PeBe ij =Pe−νe ijIm×m for all j6=i. Then, the maximum eigenvalue of Kij(uij)is non-positive, which completes the proof. Notice that it is essential to apply properly the constraints at the flux FE approximation, i.e., f0(uk) = P{i∈Nh:Cki6=0}Ckif0(ui). Otherwise, it is not possible to formally prove local bound preservation. However, experimental results in the present chapter show that using f0(uk)does not affect the overall performance of the scheme. 5.3.1 Differentiable stabilization In the case of steady, or implicit time integration, differentiability plays a role in the convergence behavior of the nonlinear solver. This is especially important if one wants to use Newton’s method. In the previous chapters we showed that nonlinear convergence can be improved after few modifications to make the scheme twice-differentiable. In this section, we introduce a set of regularizations applied to all non-differentiable functions present in the stabilized scheme introduced above. In order to regularize these functions, we follow the same strategy as in the previous chapters. Absolute values are replaced by |x|1,εh=px2+εh,|x|2,εh=x2 px2+εh . Note that |x|2,εh≤ |x| ≤ |x|1,εh. Next, we also use the smooth maximum function max σh(x, y). =|x−y|1,σh 2+x+y 2≥max(x, y). In addition, we need a smooth function to limit the value of any given quantity to one. To this end, we use Z(x). =(2x4−5x3+ 3x2+x, x < 1, 1, x ≥1.
104 Chapter 5. Monotonicity-preserving FE schemes with AMR The set of twice-differentiable functions defined above allows us to redefine the stabilization term introduced in Sect. 5.3. In particular, we define ˜ Bh(wh;uh,vh). = Pi∈NhPj∈Nh(Ωi)˜νij(wh)viuj`(i, j),for m= 1, PKe∈ThPi,j∈Nh(Ke) 1≤β,γ≤m ˜νe ij(wh)`(i, j)vβ i·δβγuγ j,for m > 1, (5.11) where ˜νij(wh). = max σh{αεh,i(wh)Kij,0, αεh,j(wh)Kji}for j∈ Nh(Ωi)\{i}, ˜νii(wh). =X j∈Nh(Ωi)\{i} ˜νij(wh), and ˜νe ij(wh). = max σhαεh,i(wh)λmax ij ,αεh,j(wh)λmax ji +X k∈M(i) Cki max σhαεh,k(wh)λmax kj ,αεh,j(wh)λmax jk +X k∈M(j) Ckj max σh(αεh,i(wh)λmax ik ,αεh,k(wh)λmax ki ) +X k∈M(i)∩M(j) CkiCkjαεh,k(wh)λmax kk ,for j∈ Nh(Ωi)\{i}, ˜νe ii(wh). =X j∈Nh(Ωi)\{i} ˜νe ij(wh). Let us note that λmax ij needs to be regularized as λmax ij =vijce ij1,εh +ckce ijk. The shock detector is also regularized as follows: αεh,i(uh). = max σh{αεh,i(uβ h)}β∈C. In the case of the component shock detector we recall the definition in Chapter 2 αεh,i(uβ h). = Z Pj∈Nh(Ωi)r∇uβ hzij1,εh +ζh Pj∈Nh(Ωi)2∇uβ h·ˆ rij2,εhij +ζh q ,(5.12) where ζhis a small value for preventing division by zero. Finally, the twice-differentiable stabilized scheme reads: find uh∈Vhsuch that uh=uhon ∂Ω,uh=u0hat t= 0, and ˜ M(un+1 h)δtUn+1 +˜ Kij(un+1 h)Un+1 =Gfor n= 1, ..., nts,(5.13)
5.4. Adaptive mesh refinement 105 where ˜ Mij(un+1 h). = [1 −max σh(αεh,i,αεh,j)] Mij + max σh(αεh,i,αεh,j)δij X j Mij, ˜ Kij(un+1 h). =Kij(un+1 h) + ˜ Bij(un+1 h), and ˜ Bij(uh) = ˜ Bh(uh;ϕj, ϕi), for i, j ∈ Nh. Corollary 5.3.3. The scheme in (5.3) with the differentiable stabilization in (5.11) is local bounds preserving, as defined in Def. 5.2.4, at any region where uβ hhas extreme values for every βin C. Proof. For an extreme value of uβ h, since |x|2,εh≤ |x|≤|x|1,εhthe quotient of (5.12) is larger than one. Hence, by definition of Z(x),αεh,i is equal to 1. At this point, it is easy to check that ˜νe ij ≥νe ij in virtue of the definition of max σh. Therefore, ρ(˜ Be ij(uh)) ≥ ρ(Be ij(uh)), completing the proof. Moreover, it is important to mention that the differentiable shock detector is weakly linearly-preserving as ζhtends to zero. This result follows directly from Sect. 2.7. In order to obtain a differentiable operator, we have added a set of regularizations that rely on different parameters, e.g., σh, εh, ζh. Giving a proper scaling of these parameters is essential to recover theoretic convergence rates. In particular, we use the following relations σh=σ|λmax|2L2(d−3)h4, εh=εL−4h2, ζh=L−1ζ, where for the scalar problem λmax is simply kvk,dis the spatial dimension of the problem, and Lis a characteristic length. 5.4 Adaptive mesh refinement The motivation of an adaptive FE method is to solve (5.10) up to a certain tolerance (or resolution) using the minimum number of DOFs. To this end, the solution error (eh=u−uh) is estimated at each element. With this information at hand, it is possible to iteratively adapt the resolution of the mesh at certain regions. This process can be divided into two parts: estimating the error at every cell, and deciding which and how many cells need to be refined or coarsened. This procedure is performed iteratively until a desired tolerance is achieved or, alternatively, a number of elements is reached. In the present chapter, we start with a rather coarse mesh and perform the following steps till reaching a stopping criterion: 1. Compute solution uh; 2. Estimate the error eh; 3. Select all cells that need to be refinement or coarsened;
106 Chapter 5. Monotonicity-preserving FE schemes with AMR 4. Update the mesh, and project the solution to the new mesh. In some cases, the refinement might be driven by features of the solution instead of a classical error estimator. For instance, one may decide to refine the regions around discontinuities. In this scenario, one could use a expression that does not estimate the error, but it allows to concentrate the elements around discontinuities. 5.4.1 Error estimators One of the keys of AMR is the ability to provide a good estimation of the error. Several error estimators have been proposed to date [1, 33, 50, 52, 94]. These can be classified, at least, in two main types. Some authors [33, 50, 78, 79, 87] try to compute an upper bound of the error for every cell. Then, provided a user defined tolerance, one can decide to refine or coarsen each cell. However, an adjoint problem needs to be solved in order to compute this upper bound [22, 50]. It is possible to approximate the error bounds without solving an adjoint problem only for simple cases, see [50]. Therefore, this kind of error estimators increases the computational cost substantially. Alternatively, one can simply determine the distribution of the error in the mesh and use this information to drive an adaptivity algorithm. In this scenario, some authors [1, 15, 63, 77, 94, 95] drive the adaptivity process with the solution gradient. In this case, explicit expressions of the estimated error are possible, requiring less computational resources than the previous option. In general, the adaptive procedure can be described as follows. Given a finite element solution uh, the error ehis approximated as eh≈∇u−∇uh. Then, the reconstruction is used as an approximation of the exact gradient. This strategy is based on superconvergence of special recovery techniques (see [95] and refs. in [1, 77]). Kuzmin and co-workers [15, 77] follow [94] to reconstruct an approximation of the exact gradient. Kelly et al. [52] proposed a well-known estimator based on gradient recovery: η2 K . =hK 24 Z∂K s∂uh ∂n {2 dΓ, where ηKis the estimated error at every element, K. The main advantages of this estimator are its simplicity and its low computational cost. For these reasons, this estimator is used in the present chapter. It is worth mentioning that our problems of interest are characterized by exhibiting discontinuities, where the error concentrates. These regions are susceptible to develop instabilities, and thus, these are the regions in which the shock capturing is activated. Therefore, it is natural to use the shock capturing to drive the adaptivity procedure. We propose an estimator based on the graph Laplacian `(i, j)present in the stabilization term (5.4). This way, we reuse available information and reduce the computational
5.5. Nonlinear solver 107 overhead associated with error estimation. The indicator reads: ˜η2 K . =hd−2 K`(uβ h,uβ h)K=hd−2 KX i∈Nh(K)X j∈Nh(Ωi) (uβ i−uβ j)2, where β∈Cis the index of the specific component analyzed. This estimator is expected to yield high values around shocks and low values in smooth regions. 5.4.2 Refinement strategy After the error has been estimated for every element, one needs to decide which element needs to be refined and which one coarsened. If an upper bound of the error is computed, then one may use a given tolerance to make this decision. However, in the present case this is not available. A classical alternative is to refine/coarsen a fixed amount of elements at every iteration [10, 12]. In the present study, a 30% of the elements with higher error estimates are refined whereas a 10% of the elements with lower error estimates are coarsened. This percentages are arbitrary and other choices are valid. Notice that using this setting in two dimensions the number of elements is almost doubled at every iteration. We make use of the parallel nth element algorithm [10, 90] to efficiently determine the error estimator thresholds for refining or coarsening the elements. 5.5 Nonlinear solver In this section, we describe the method used for solving the nonlinear system of equations arising from the scheme introduced above. In particular, we use a hybrid Picard–Newton approach in order to increase the robustness of the nonlinear solver. Moreover, we also make use of a line-search method to improve the nonlinear convergence. We define the residual of the equation (5.13) at the k-th iteration as R(uk,n+1 h). =˜ M(uk,n+1 h)δtUk,n+1 +˜ Kij(uk,n+1 h)Uk,n+1 −G. Hence, the Jacobian is defined as J(uk,n+1 h). =∂R(uk,n+1 h) ∂Uk,n+1 (5.14) =˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h) + ∂˜ M(uk,n+1 h) ∂Uk,n+1 δtUk,n+1 +∂˜ Kij(uk,n+1 h) ∂Uk,n+1 Uk,n+1. Therefore, Newton method consists in solving J(uk,n+1 h)∆Uk+1,n+1 =−R(uk,n+1 h). It is well known that Newton method can diverge if the initial guess of the solution u0,n+1 h is not close enough to the solution. In order to improve robustness, we use a linesearch method to update the solution at every time step. The new approximation is computed as Uk+1,n+1 =Uk,n+1 +λ∆Uk+1,n+1, where λis obtained using a standard cubic backtracking algorithm.
108 Chapter 5. Monotonicity-preserving FE schemes with AMR As introduced at the beginning of the section, we also use a hybrid approach combining Newton method with Picard linearization. Picard nonlinear iterator can be obtained removing the last two terms of (5.14), i.e., ˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h)∆Uk+1,n+1 =−R(uk,n+1 h).(5.15) Clearly, it is equivalent to ˜ M(uk,n+1 h) + ˜ Kij(uk,n+1 h)Uk+1,n+1 =˜ M(uk,n+1 h)Un+G. Moreover, we modify the left hand side terms in (5.15); we use αi= 1 for computing these terms while we use the value obtained from (5.8) for the residual. Using this strategy, the solution remains unaltered but the obtained approximations uk,n+1 hfor intermediate values of kare more diffusive. Even though this modification slows the nonlinear convergence, it is essential at the first iterations. Otherwise, the robustness of the method might be jeopardized. The resulting iterative nonlinear solver consists in the following steps. We iterate Picard method in (5.15), with the modification described above, until the L2norm of the residual is smaller than a given tolerance. In the present chapter, we use a tolerance of 10−2. Afterwards, Newton method with the exact Jacobian in (5.14) is used until the desired nonlinear convergence criteria is satisfied. We summarize the nonlinear solver introduced above in Alg. 4. Algorithm 4: Hybrid Picard–Newton method. Input:U0,n+1,tol1,tol2,ε Output:Uk,n+1,k k= 1,ε1=ε while kR(Uk,n+1)k/kR(U0,n+1)k ≥ tol1do Compute αi(Uk,n+1)using (5.8) Compute ∆Uk+1,n+1 using (5.15) Minimize kR(Uk+1,n+1)k, where Uk+1,n+1 =λ∆Uk+1,n+1 +Uk,n+1, with respect to λ Set Uk+1,n+1 =λ∆Uk+1,n+1 +Uk,n+1 Update k=k+ 1 while kR(Uk,n+1)k/kR(U0,n+1)k ≥ tol2do Compute αi(Uk,n+1)using (5.8) Solve J(Uk,n+1)∆Uk+1,n+1 =−R(Uk,n+1)with Jin (5.14) Minimize kR(Uk+1,n+1)k, where Uk+1,n+1 =λ∆Uk,n+1 +Uk,n+1, with respect to λ Set Uk+1,n+1 =λ∆Uk,n+1 +Uk,n+1 Update k=k+ 1
5.6. Numerical results 109 5.6 Numerical results In this section, we perform several numerical experiments to assess the numerical scheme introduced in the previous sections. First, we perform a convergence analysis to assess its implementation. Then, we use steady benchmark tests to analyze the effectiveness of the high-order scheme in the context of AMR. In particular, we compare the nonlinear scheme in (5.13) with its linear (first order) counterpart, i.e., using αεh,i(uh)≡1. From the experience in the numerical experiments of the previous chapters, we choose the following regularization parameters: σ= 10−2,ε= 10−4, and γ= 10−10. In addition, for all Euler tests below, the density is discontinuous at all shocks. Therefore, we use C={1}in (5.8), i.e., the shock detector is based on the density behavior. 5.6.1 Convergence First, the convergence to a discontinuous solution is analyzed. To this end, we solve two different problems. On the one hand, the following scalar problem is solved ∇·(vu) = 0 in Ω = [0,1] ×[0,1], u=uDon Γin,(5.16) where v(x, y). = (1 /2,sin −π/3), and inflow boundary conditions uD= 1 on {x= 0}∩{y > 0.7}and y= 1, while uD= 0 at the rest of the inflow boundary. This problem has the following analytical solution u(x, y) = (1if y > 0.7+2xsin −π/3, 0otherwise. For the Euler equations, the problem is the well known compression corner test [3, 61], also known as oblique shock test [84, 88]. This benchmark consists in a supersonic flow impinging to a wall at an angle. We use a [0,1]2domain with a M= 2 flow at 10◦ with respect to the wall. This leads to two flow regions separated by an oblique shock at 29.3◦, see Fig. 5.3. Figure 5.3: Compression corner scheme. Since the solution is not smooth, we expect linear convergence rates in the L1-norm. Fig. 5.4 shows the convergence behavior of both problems with uniform mesh refinements.
110 Chapter 5. Monotonicity-preserving FE schemes with AMR 16 32 64 128 1/h 10−2 10−1 kuh−ukL1(Ω) 11 (a) Scalar transport problem. 16 32 64 128 1/h 3·10−3 10−2 5·10−2 kuh−ukL1(Ω) 11 (b) Euler equations. Figure 5.4: Convergence of ku−uhkL1(Ω) to a solution with a discontinuity. The experimental convergence rate measured for the scalar transport problem is 0.82, whereas the convergence rate measured is 0.94 for the compression corner test. Therefore, both tests exhibit the expected convergence behavior. 5.6.2 Linear discontinuity For this test, we use again the problem in (5.16). The purpose of this test is twofold. On the one hand, we analyze the effectiveness of the proposed estimator. On the other hand, we compare the effectiveness of the linear and nonlinear stabilization methods. Specifically, this effectiveness is measured as follows. For a given error, we consider a method more effective if it requires less computational time, independently of the number of elements required. In addition, we also solve the problem for successive uniformly refined meshes in order to evaluate the effect of AMR. For all comparisons, we start with a coarse mesh of 16 ×16 elements, and proceed adapting the mesh up to a maximum number of elements. For the nonlinear stabilization, we set a maximum of 104elements. The maximum number of elements for the low-order method is 105. The uniform mesh is refined up to a 256 ×256 mesh. We use a nonlinear tolerance of k∆uhk/kuhk<10−4, and a maximum of 500 iterations. Fig. 5.5 shows the evolution of the AMR algorithm for both estimators. The results shown in this picture have been obtained using the linear stabilization, and the leftmost column using the nonlinear one. It can be observed that both Kelly (ηK) and graph Laplacian (˜ηK) estimators refine in the vicinity of the shock. However, the graph Laplacian operator clearly outperforms Kelly estimator. Figs. 5.6–5.8 compare the effectiveness of the low-order and the high-order stabilization schemes. The results are obtained for the stabilization parameter q= 1,q= 2, and q= 10, respectively.
5.6. Numerical results 117 In Fig. 5.13, we depict the refinement evolution for the graph Laplacian estimator (˜ηK) for linear and nonlinear stabilization. As expected, we can observe that for the high-order method the scheme is able to resolve the shock with less refinement steps. The linear stabilization is able to provide well-resolved shocks at the final refinement step. Fig. 5.14 compares the effectiveness of the low-order and the high-order stabilization schemes for different values of q. The high-order scheme is able to converge efficiently and the overhead of solving a nonlinear problem does not affect the overall performance. In this case, the low-order and the high-order schemes require similar computational time for any given error. Actually, for the finer meshes, the high-order scheme with either q= 1 or q= 2 already performs better than the low-order scheme. However, for some meshes the high-order scheme exhibits convergence problems. In the case of q= 10, the cost of converging the nonlinear problem does not compensate the increase in computational cost. 10−1100101102 Time [s] 10−2 10−1 kρh−ρkL1(Ω) 103104105 Number elements 10−2 10−1 kρh−ρkL1(Ω) Low order, ˜η High order, ˜η,q= 1 High order, ˜η,q= 2 High order, ˜η,q= 10 Figure 5.14: Time and elements convergence comparison for the compression corner problem. 5.6.5 Reflected shock This benchmark consists in two flow streams colliding at different angles. The domain has dimensions [0.0,1.0]×[0.0,4.1] and a solid wall at its lower boundary. This configuration leads to a steady shock separating both flow regimes that is reflected at the wall producing a third different flow state behind it. A sketch of this benchmark test is given in Fig. 5.15. The flow states at each region have been collected in Tab. 5.1. Table 5.1: Reflected shock solution values at every region. Region Density [Kg m−3] Velocity [m s−1] Total energy [J] a 1.0 (2.9, 0.0) 5.99075 b 1.7 (2.62, -0.506) 5.8046 c 2.687 (2.401, 0.0) 5.6122
118 Chapter 5. Monotonicity-preserving FE schemes with AMR Figure 5.15: Reflected shock scheme. We analyze the effectiveness of the high-order scheme, and evaluate the performance of the graph Laplacian estimator. We start with a coarse mesh of 16 ×64 elements and adapt the mesh till a certain number of elements is reached. For the high-order method, we set a maximum of 104elements. The maximum number of elements for the low-order method is 3·105. We use a nonlinear tolerance of k∆uhk/kuhk<10−4and a maximum of 500 iterations. Fig. 5.16 compares the effectiveness of the low-order and the high-order stabilization schemes for different values of q. The high-order scheme converges efficiently and the overhead of solving a nonlinear problem does not affect the overall performance. Actually, for the most refined meshes the high-order method is more efficient than the low-order one. As for the previous problem, Fig. 5.16 shows that the high-order scheme can present nonlinear convergence problems at some steps of the refinement process. However, as the mesh becomes more adapted to the problem this issues is reduced. 100101102103 Time [s] 10−1 100 kρh−ρkL1(Ω) 103104105 Number elements 10−1 100 kρh−ρkL1(Ω) Low order, ˜η High order, ˜η,q= 1 High order, ˜η,q= 2 High order, ˜η,q= 10 Figure 5.16: Time and elements convergence comparison for the reflected shock problem. In Fig. 5.17 we depict the refinement evolution for the graph Laplacian estimator (˜ηK) for the low-order scheme. In these figures it can be observed how the graph Laplacian estimator is able to concentrate all the resolution at the shock location. Finally, we can conclude from the lower two figures that both schemes resolve the shocks properly after the mesh has been refined enough.
5.6. Numerical results 119 Figure 5.17: Evolution of the mesh refinement process. ˜ηKwith loworder scheme is used. For the low-order scheme from top to bottom results have been obtained at refinement step 1, 2, 3, 4, 5, 6, and 7. The lower two figures are the high-order with q= 2 (top) and low-order (bottom) results at their last refinement step.
120 Chapter 5. Monotonicity-preserving FE schemes with AMR 5.7 Conclusions The stabilization schemes in Chapters 2 and 4 have been extended and assessed in the AMR context for nonconforming hierarchical octree meshes. The chapter focuses in assessing the effectiveness of linear (first-order) and nonlinear (higher-order) stabilization. We focus the comparison in terms of accuracy versus computational time. The results indicate that linear stabilization is more effective for coarse meshes. In this case, the computational cost required to solve the stiff nonlinear problem due to the nonlinear stabilization does not compensate the improvement in the accuracy. This is especially evident for linear systems of PDEs. On the contrary, as the mesh is refined and properly adapted to the shocks, nonlinear stabilization pays the price. Even though increasing the value of qin the nonlinear stabilization (a parameter that makes shocks sharper but hinders nonlinear convergence) improves accuracy, it turns to be far more effective to refine the mesh further for low values of q. Nevertheless, it is worth mentioning that high-order method might exhibit nonlinear convergence problems for some meshes. In addition, a new refinement criterion have been proposed. The proposed estimator is based on the graph Laplacian used in the definition of the stabilization method. Numerical results show that this shock detector is able to perform better that the well known Kelly estimator for problems with shocks or discontinuities.
Chapter 6 Conclusions and future work 6.1 Conclusions In this thesis, the development of monotonicity-preserving FE methods has been explored. Since the main chapters of this dissertation are self-contained and preserve the structure of a paper, each one contains their own detailed conclusions. In this chapter, we present a more general overview. To this end, let us recall the list of goals set in Sect. 1.2. •Design of a monotonicity-preserving scheme for arbitrary mesh geometries. In Chapter 2 we consider a nonlinear stabilization technique for the FE approximation of scalar conservation laws with implicit time stepping. The method relies on an artificial diffusion method, based on a graph-Laplacian operator. The artificial diffusion term is based on a graph-Laplacian artificial diffusion operator, instead of a PDE-based one. This removes any requirements on the mesh. In addition, this strategy is also used in subsequent chapters. The stabilization method proposed in Chapter 2 satisfies the local DMP, and thus preserves monotonicity. Furthermore, all numerical results in Sect. 2.9 exhibit this property. •Analysis and improvement of the nonlinear convergence behavior of monotonicity-preserving schemes. The scheme proposed in Chapter 2 is proved to be Lipschitz continuous. This property leads to well-posedness of the nonlinear problem. However, the resulting scheme is highly nonlinear, leading to very poor nonlinear convergence rates. Therefore, we also propose a regularized version of the scheme that is twice differentiable. This allowed us to use Newton’s method with the exact Jacobian. Numerical experiments in Sect. 2.9 show a reduction of 10 to 20 times in the number of iterations with respect to the original non-differentiable algorithms. •Extension to high-order discretizations in space and time. In Chapter 3, the stabilization method in Chapter 2 is extended to isogeometric analysis. The proposed method is DMP-preserving for arbitrary high-order discretizations in space and time without any CFL-like condition. Moreover, in 121
122 Chapter 6. Conclusions and future work order to reduce the computational cost of the space–time method, we propose a partitioned scheme in Sect. 3.4. This alternative scheme is also proved to be unconditionally DMP-preserving. In addition, all numerical in results Sect. 3.6 exhibit this property. •Extension to first order hyperbolic systems of equations. Chapter 4 is devoted to this goal. Extension of monotonicity preservation to systems of equations is not straightforward. Actually, the continuous problem does not necessarily need to have monotonic solutions. In this case, previous works in literature resort to proving local bounds preservation. This could be seen as an heuristic extension of the properties required to prove the LED property for scalar problems (see Th. 2.4.1 and Sect. 4.2.3). In Chapter 4, a differentiable local bounds preserving stabilization method for Euler equations is developed. In addition, a continuation method for the regularization parameters present in the differentiable stabilization is proposed to improve nonlinear convergence. Numerical results show that differentiability not only can improve nonlinear convergence, but it also improves the robustness of the method. However, the improvement in the nonlinear convergence is restricted to moderate tolerances, and is not as significant as in Chapter 2. •Extension to AMR FE schemes. In Chapter 5, the stabilization schemes in Chapter 2 and Chapter 4 are extended and assessed in the AMR context. In particular, we use nonconforming hierarchical octree meshes. The stabilization method for scalar problems is defined using the assembled matrix. Thus, it can work directly with this kind of meshes without any modification. Instead, the stabilization method for systems of equations uses elemental values to define the artificial diffusion. Hence, minor modifications are introduced to adapt the scheme. In any case, with these modifications, the scheme preserves all the properties of the original method (see 5.3.3). Moreover, a new refinement criterion is proposed in Sect. 5.4. The proposed estimator is based on the graph Laplacian used in the definition of the stabilization method. Numerical results show that this shock detector is able to perform better than the well known Kelly estimator for problems with shocks or discontinuities. •Assessment of the efficiency of high-order monotonicity-preserving schemes in AMR context. The results of Chapter 5 indicate that the high-order scheme only becomes superior for a sufficiently refined mesh. The low-order method might perform better for coarse meshes. However, the high-order scheme has a higher convergence rate. Thus, it outperforms the low-order scheme once the mesh is refined enough.
6.2. Future work 123 In addition, it can be observed in Sect. 5.6 that it is more efficient to refine the mesh than to improve accuracy by using high values of q. Let us recall that qis a parameter in the stabilization that allows one to modulate the amount of artificial diffusion introduced. The higher it is, the less diffusive the presented method is. Nevertheless, it is worth mentioning that the high-order method might exhibit nonlinear convergence problems for some meshes. In this scenario, the low-order scheme clearly outperforms the high-order method. 6.2 Future work Research never comes to an end, it is simply bounded by time. In this section, we proceed to describe a few ideas that arise as possible continuation of the developments in this dissertation. •Smoothness indicator for isogeometric analysis The high-order method developed in Chapter 3 yields solutions that satisfy the global DMP for arbitrary order discretizations. However, only for monotonic solutions the method is able to recover the high-order convergence rates. This is a direct consequence of the stabilization method introduced, which is only second order accurate. Lohmann et al. [71] proposed to use smoothness indicators based on the behavior of second order derivatives. This prevents to formally proof monotonicity preservation, but numerical experiments show an improved behavior and high-order convergence rates can be recovered. An interesting work could be to develop such a smoothness indicator for the particular case of isogeometric analysis. In this case, one could take advantage of the higher continuity of this kind of discretizations. •Parallelization of the implementation The methods presented in this dissertation have only been tested in serial experiments. In any case, all methods are local and nothing prevents its parallelization. However, the domain of dependence is slightly larger than the one for regular FE methods. Therefore, the parallel implementation should be adapted to support at least one layer of ghost elements at subdomain interfaces. This is especially important if one is willing to compute the exact Jacobian. Otherwise, if an inexact Jacobian is sufficient, as performed in [18], then the parallel implementation does not require any special requirement. •AMR for high-order methods In Chapter 5, we have adapted the methods in Chapter 2 and 4 to adaptive meshes. However, this is not the case for the method in Chapter 3 due to the higher coupling of isogeometric analysis basis functions. Bornemann and Cirak [20] use hierarchical B-splines to achieve an AMR discretization. We consider interesting to
124 Chapter 6. Conclusions and future work combine this kind of B-spline discretizations with the AMR methods in Chapter 5. Furthermore, developing an hp-adaptive method using this strategy could lead to an improved behavior for problems that combine shocks with regions where the solution is smooth. •Extension to other problems The most complex problem solved in this thesis are the Euler equations. As motivated in Chapter 1, we would like to eventually extend these methods to enhance plasma simulations. Therefore, the immediate development to be performed is the extensions to compressible Navier-Stokes, and ideal magnetohydrodynamics (MHD). In a latter stage, extensions to resistive MHD, and multi-fluid plasma equations should also be performed. •Extension to compatible discretizations In combination with the previous point, we consider interesting to explore extensions to compatible discretizations. The magnetic field in MHD formulations is solenoidal. Several strategies have been developed to deal with this constraint. A common approach is to use Nédélec FEs [80] to discretize the magnetic field. Therefore, we consider that extending the methods presented in this dissertation to compatible discretizations could increase their applicability.
Bibliography [1] M. Ainsworth and J. Tinsley Oden,A posteriori error estimation in finite element analysis, Computer Methods in Applied Mechanics and Engineering, 142 (1997), pp. 1–88. [2] R. Anderson, V. Dobrev, T. Kolev, D. Kuzmin, M. Quezada de Luna, R. Rieben, and V. Tomov,High-order local maximum principle preserving (MPP) discontinuous Galerkin finite element method for the transport equation, Journal of Computational Physics, 334 (2017), pp. 102–124. [3] J. D. Anderson Jr.,Modern Compressible Flow, McGraw-Hill, 2nd ed., 1990. [4] S. Badia and J. Bonilla,Monotonicity-preserving finite element schemes based on differentiable nonlinear stabilization, Computer Methods in Applied Mechanics and Engineering, 313 (2017), pp. 133–158. [5] S. Badia, J. Bonilla, and A. Hierro,Differentiable monotonicity-preserving schemes for discontinuous Galerkin methods on arbitrary meshes, Computer Methods in Applied Mechanics and Engineering, 320 (2017), pp. 582–605. [6] S. Badia, J. Bonilla, S. Mabuza, and J. N. Shadid,Differentiable local bounds preserving stabilization for first order hyperbolic problems, Submitted, (2019). [7] S. Badia and A. Hierro,On Monotonicity-Preserving Stabilized Finite Element Approximations of Transport Problems, SIAM Journal on Scientific Computing, 36 (2014), pp. A2673–A2697. [8] S. Badia and A. Hierro,On discrete maximum principles for discontinuous Galerkin methods, Computer Methods in Applied Mechanics and Engineering, 286 (2015), pp. 107–122. [9] S. Badia and A. F. Martín,A tutorial-driven introduction to the parallel finite element library FEMPAR v1.0.0, (2019). [10] S. Badia, A. F. Martín, E. Neiva, and F. Verdugo,A generic finite element framework on parallel tree-based adaptive meshes, Submitted, (2019). [11] S. Badia, A. F. Martín, and J. Principe,FEMPAR: An Object-Oriented Parallel Finite Element Framework, Archives of Computational Methods in Engineering, 25 (2018), pp. 195–271. 125
126 BIBLIOGRAPHY [12] W. Bangerth, C. Burstedde, T. Heister, and M. Kronbichler,Algorithms and data structures for massively parallel generic adaptive finite element codes, ACM Trans. Math. Softw., 38 (2012), pp. 14:1–14:28. [13] G. R. Barrenechea, E. Burman, and F. Karakatsani,Edge-based nonlinear diffusion for finite element approximations of convection-diffusion equations and its relation to algebraic flux-correction schemes, Numerische Mathematik, (2016), pp. 1–25. [14] G. R. Barrenechea, V. John, and P. Knobloch,Analysis of Algebraic Flux Correction Schemes, SIAM Journal on Numerical Analysis, 54 (2016), pp. 2427– 2451. [15] M. Bittl and D. Kuzmin,An hp-adaptive flux-corrected transport algorithm for continuous finite elements, Computing, 95 (2013), pp. 27–48. [16] J. Bonilla and S. Badia,Maximum-principle preserving space-time isogeometric analysis, Computer Methods in Applied Mechanics and Engineering, 354 (2019), pp. 422–440. [17] , Monotonicity-preserving finite element schemes with adaptive mesh refinement for hyperbolic problems, In preparation, (2019). [18] J. Bonilla, S. Mabuza, J. N. Shadid, and S. Badia,On Differentiable Linearity and Local Bounds Preserving Stabilization Methods for First Order Conservation Law Systems, in Center for Computing Research Summer Proceedings 2018, A. Cangi and M. L. Parks, eds., Sandia National Laboratories, 2018, pp. 107–119. [19] P. Bonoli and L. C. McInnes,Report of the Workshop on Integrated Simulations for Magnetic Fusion Energy Sciences, tech. report, 2015. [20] P. B. Bornemann and F. Cirak,A subdivision-based implementation of the hierarchical b-spline finite element method, Computer Methods in Applied Mechanics and Engineering, 253 (2013), pp. 584–598. [21] S. C. Brenner and L. R. Scott,The Mathematical Theory of Finite Element Methods, vol. 15 of Texts in Applied Mathematics, Springer New York, New York, NY, softcover ed., nov 2008. [22] E. Burman,Adaptive finite element methods for compressible flow, Computer Methods in Applied Mechanics and Engineering, 190 (2000), pp. 1137–1162. [23] , On nonlinear artificial viscosity, discrete maximum principle and hyperbolic conservation laws, BIT Numerical Mathematics, 47 (2007), pp. 715–733. [24] , A monotonicity preserving, nonlinear, finite element upwind method for the transport equation, Applied Mathematics Letters, 49 (2015), pp. 141–146.