Advances in the Discontinuous Galerkin Method: Hybrid Schemes and Applications to the Reactive Infiltration Instability in an Upwelling Compacting Mantle by Alan R. Schiemenz B.S., Baylor University, Texas, 2004 Sc.M., Brown University, Rhode Island, 2005 A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the Division of Applied Mathematics at Brown University PROVIDENCE, RHODE ISLAND May 2009 Abstract of “Advances in the Discontinuous Galerkin Method: Hybrid Schemes and Applications to the Reactive Infiltration Instability in an Upwelling Compacting Mantle” by Alan R. Schiemenz, Brown University, May 2009 Abstract High-order methods are emerging in the scientific computing community as superior alternatives to the classical finite difference, finite volume, and continuous finite element methods. The discontinuous Galerkin (DG) method in particular combines many of the positive features of all of these methods. This thesis presents two projects involving the DG method. First, a Hybrid scheme is presented, which implements DG areas where the so- lution is considered smooth, while dropping the order of the scheme elsewhere and implementing a finite volume scheme with high-order, non-oscillatory solution recon- structions suitable for unstructured mesh. Two such reconstructions from the ENO class are considered in the Hybrid. Successful numerical results are presented for nonlinear systems of conservation laws in one dimension. Second, the high-order discontinuous Galerkin and Fourier spectral methods are applied to an application modeling three-phase fluid flow through a porous medium, undergoing solid-fluid reaction due to the reactive infiltration instability (RII). This model incorporates a solid upwelling term and an equation to track the abundance of the reacting mineral orthopyroxene (opx). After validating the numerical discretiza- tion, results are given that provide new insight into the formation of melt channels in the Earth’s mantle. Mantle heterogeneities are observed to be one catalyst for the development of melt channels, and the dissolution of opx produces interesting bifurcations in the melt channels. An alternative formulation is considered where the mass transfer rate relative to velocity is taken to be infinitely large. In this setting, the stiffest terms are removed, greatly reducing the cost of time integration. c Copyright 2009 by Alan R. Schiemenz This dissertation by Alan R. Schiemenz is accepted in its present form by the Division of Applied Mathematics as satisfying the dissertation requirement for the degree of Doctor of Philosophy. Date Jan S. Hesthaven, Ph.D., Advisor Recommended to the Graduate Council Date Chi-Wang Shu, Ph.D., Reader Date Edgar M. Parmentier, Ph.D., Reader Approved by the Graduate Council Date Sheila Bonde, Dean of the Graduate School iii The Vita of Alan R. Schiemenz Born on August 12, 1983 in Waco, Texas. Education • Sc. M. in Applied Mathematics, Brown University, May 2005. • B.S. in Mathematics, Baylor University, Waco, Texas, May 2004. Publications/Preprints • Schiemenz, A., Hesse, M., Liang, Y., Parmentier, E., and Hesthaven, J., Reac- tive Infiltration Instability in an Upwelling Compacting Mantle: Linear Stabil- ity Analysis, in progress. • Schiemenz, A., Hesse, M., Liang, Y., Parmentier, E., and Hesthaven, J., Re- active Infiltration Instability in an Upwelling Compacting Mantle: Nonlinear Numerical Simulations, in progress. • Schiemenz, A. and Robinson, A., CSRI Summer Proceedings 2007 - The Com- puter Science Research Institute at Sandia National Laboratories, edited by M. L. Parks and S. S. Collis, http://www.cs.sandia.gov/CSRI/Proceedings/CSRI2007.pdf, SAND2007-7977, 2007, pp. 231-241, Energy Based Magnetic Forces Compu- tation using Automatic Differentiation. Research Interests • Numerical solutions to partial differential equations. • Finite element methods, particularly the discontinuous Galerkin (DG) method. • High-order numerical methods for fluid flow through porous media. • High-order limiting and essentially non-oscillatory (ENO) reconstruction meth- ods. iv Academic Experience • Teaching Assistant: Methods of Applied Mathematics, I, Division of Applied Mathematics, Brown University, Fall 2005 & 2008. • Grading Assistant: Numerical Solutions to Partial Differential Equations, II, Division of Applied Mathematics, Brown University, Spring 2008. • Private Mathematics Tutor, January 2007 - May 2009. • Math Resource Center, Brown University, August 2006 - May 2007. Awards and Honors • Student Travel Support Grant, SIAM conference on mathematical & compu- tational issues in the geosciences, June 2009, Leipzig, Germany. • Dissertation Fellowship, Brown University, Spring 2009. Graduate Student Fellowship, Division of Applied Mathematics, Brown University, 2003-2004. • Phi Beta Kappa, inducted 2004 into the Zeta of Texas chapter. • Outstanding Mathematics Student, 2004, Baylor University. Presentations • SIAM Conference on Mathematical & Computational Issues in the Geosciences, June, 2009. High-order Methods for the Reactive Infiltration Instability in an Upwelling Compacting Mantle. • American Geophysical Union Fall Meeting, December, 2008. A Possible Link between Mantle Heterogeneity and Dunite Channel Distribution in an Up- welling Mantle: A Numerical Study. • Sandia Student Symposium, August 2007. Improving Magnetic Energy Force Discretization in ALEGRA. v • Brown University, Division of Applied Mathematics “Pizza Seminar” Series,May 2007. Discontinuous Galerkin Method: Introduction, Analysis and use in Hy- brid Schemes. vi Acknowledgments Foremost I would like to thank the Lord for His providential plan in my life. When I arrived at Brown in the fall of 2004 I set foot in a new world, with little understanding of what awaited. In the five years since much has been revealed, and through the enduring experience of graduate school the Lord has granted me the opportunity for growth both in my professional endeavors and my character. I would like to sincerely thank my advisor, Professor Jan S. Hesthaven, for his willingness to take me as a student. He has invested a great amount of his time and effort in my education, always providing friendly mentorship and encouragement. My thanks go out to my colleagues, Dr. Marc Hesse, and Professors Jan Liang and Marc Parmentier. They have introduced me to the field of computational geo- sciences. Moving to the application realm of numerical analysis has greatly shaped my professional interests. I would like to express my appreciation to my friends and fellow graduate stu- dents, Akil Naryan, Andreas Kloeckner, Michael Lamar, Ravi Srinivasan and Kenny Chowdhary. They have provided me with great friendship and technical assistance over the years. Finally, I am forever indebted to my family. I thank my parents for their unwa- vering support in often-difficult times during my graduate career. I thank my wife for her grace and companionship during this time in my life. Without the encour- agement and love of my family this accomplishment would have been impossible. vii Contents Abstract Acknowledgments vii 1 Discontinuous Galerkin Methods 1 1.1 History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Basic Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Explicit Runge-Kutta Time Integration . . . . . . . . . . . . . . . . . 6 1.5 Problems With Nonlinear Problems . . . . . . . . . . . . . . . . . . . 8 1.6 Introduction to Finite Volume Schemes . . . . . . . . . . . . . . . . . 9 1.6.1 Numerical Flux . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 I Hybrid Discontinuous Galerkin Methods 13 2 ENO-Type Reconstructions 14 2.1 Properties of the ENO Class . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.1 Smoothness Indicators . . . . . . . . . . . . . . . . . . . . . . 17 2.2 QENO Reconstruction Scheme . . . . . . . . . . . . . . . . . . . . . . 18 2.2.1 Conservation of Mean . . . . . . . . . . . . . . . . . . . . . . 18 2.2.2 Stencil Selection . . . . . . . . . . . . . . . . . . . . . . . . . . 18 viii 2.2.3 Reconstruction Algorithm . . . . . . . . . . . . . . . . . . . . 19 2.2.4 Weights and the ENO Property . . . . . . . . . . . . . . . . . 21 2.3 DK Reconstruction Scheme . . . . . . . . . . . . . . . . . . . . . . . 22 2.3.1 Polynomial Representation . . . . . . . . . . . . . . . . . . . . 22 2.3.2 Stencil Selection . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.3 Reconstruction Algorithm . . . . . . . . . . . . . . . . . . . . 24 2.3.4 Weights and the ENO Property . . . . . . . . . . . . . . . . . 28 3 The Hybrid Scheme 31 3.1 Switching Between Schemes . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.1 Implementation Algorithm . . . . . . . . . . . . . . . . . . . . 33 3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.1 Convergence Studies . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.2 One-dimensional Advection Equation . . . . . . . . . . . . . . 36 3.2.3 Burgers’ Equation . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.2.4 One-dimensional Euler Equations . . . . . . . . . . . . . . . . 46 3.3 Two-Dimensional Observations . . . . . . . . . . . . . . . . . . . . . 53 3.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 II Reactive Infiltration Instability 58 4 High-Order Methods for the Reactive Infiltration Instability in an Upwelling Compacting Mantle 59 4.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.1.1 Reaction Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.1.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Dimensionless Equilibrium Equations on a Reactive Time Scale . . . 68 4.3 One-Dimensional Steady State . . . . . . . . . . . . . . . . . . . . . . 70 4.3.1 Dimensionless Boundary Conditions . . . . . . . . . . . . . . . 72 ix 4.3.2 Domain Length . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4 Numerical Implementation . . . . . . . . . . . . . . . . . . . . . . . . 73 4.4.1 Fourier Spectral Discretization . . . . . . . . . . . . . . . . . . 75 4.4.2 Discontinuous Galerkin Discretization . . . . . . . . . . . . . . 76 4.4.3 Representing the Solution . . . . . . . . . . . . . . . . . . . . 82 4.4.4 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . 83 4.4.5 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . 87 4.4.6 Additional Remarks . . . . . . . . . . . . . . . . . . . . . . . 87 5 Validation of Numerical Model 89 5.1 Linear Stability Analysis for Kinetic Formulation . . . . . . . . . . . 89 5.2 Linear Stability Analysis for Equilibrium Formulation . . . . . . . . . 91 5.3 Convergence of Kinetic to Equilibrium Formulations . . . . . . . . . . 92 5.4 Matching Transient Growth Rate to Linear Stability Analysis . . . . 93 6 Observations from Numerical Model 96 6.1 Effects of the φm f Viscosity Term . . . . . . . . . . . . . . . . . . . . . 96 6.2 Effects of the φlopx Surface Area Term . . . . . . . . . . . . . . . . . . 97 6.3 Dimensionless Governing Equations on an Upwelling Time Scale . . . 98 6.4 Transient Numerical Results for the Kinetic Formulation . . . . . . . 101 6.4.1 Correlation of Inflow φf - Perturbations to Outflow Melt Channels101 6.4.2 Channels Deepen as Da → ∞ . . . . . . . . . . . . . . . . . . 104 6.4.3 Observations From a Single Melt Channel . . . . . . . . . . . 105 6.5 Transient Numerical Results for the Equilibrium Formulation . . . . . 107 6.6 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 x List of Tables 3.1 Convergence study for DK and Q-ENO reconstructions. L∞ error is given for K elements of order N . . . . . . . . . . . . . . . . . . . . . 35 3.2 Convergence study for DK and Q-ENO reconstructions. L1 error is given for K elements of order N . For DK we use a biasing parameter of λ1 = 103 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3 Convergence study of Burgers’ equation with smooth initial condition u(x, 0) = − sin (πx). L∞ errors are computed at time t = 0.2. Results are given using a A pure finite volume scheme with Q-ENO and DK reconstructions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.4 Convergence study of Burgers’ equation with smooth initial condition u(x, 0) = − sin (πx). L1 errors are computed at time t = 0.2. A pure finite volume scheme with Q-ENO and DK reconstructions is used . . 43 3.5 Two-dimensional convergence study for DK and Q-ENO reconstruc- tions for the function u(x, y) = exp [−5 (x2 + y 2 )]. L1 error is given for K elements of order N . The schemes attain their desired recon- struction order. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Typical values for the dimensionless parameters given in (4.32). . . . 69 6.1 Typical values for the dimensionless parameters given in (6.2). . . . . 99 xi List of Figures 3.1 One-dimensional advection, Riemann problem (see Eq. (3.8)). Termi- nal time is t = 1. O’s: Q-ENO reconstruction. X’s: DK reconstruc- tion. Solid line: exact solution. Left to right: K = 20, 40, 80. Top to bottom: N = 1, 2, 3, 4. The order of polynomial reconstruction is the same as the order of the DG scheme. A switch tolerance of 0 = 10−4 is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2 Hybrid solution for one-dimensional advection equation with initial condition given by Eq. (3.9). Q-ENO reconstruction is given by O’s; DK is given by X’s. Exact solution is given by solid line. K = 200 elements are used, with varying orders N = 1, 2, 3, 4. Switch tolerance is fixed at 0 = 10−4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 3.3 Hybrid solution for one-dimensional advection equation with initial condition given by Eq. (3.9). Q-ENO reconstruction is given by O’s; DK is given by X’s. Exact solution is given by solid line. A mesh of K = 200 elements are used, with varying order N = 1, 2, 3, 4. Switch tolerance is 0 = 10−2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.4 Solution to Burgers’ equation (3.11)-(3.12) at times t = 0.2 (left) and t = 1 (right). A uniform mesh of order N = 3 and K = 80 elements is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.5 Error plots corresponding to the simulations of Burgers’ equation in Figure 3.4 are given. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 xii 3.6 Solution to Burgers’ equation (3.11)-(3.12) at times t = 0.2 (left) and t = 1 (right). A non-uniform mesh of order N = 4 and K = 40 elements is used, with element stretching corresponding to Eq. (3.15). 45 3.7 Numerical solution to Sod’s problem (3.18). K = 200 elements and N = 1 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted. O’s: Q-ENO reconstruction. X’s: DK recon- struction. Solid line: exact solution. . . . . . . . . . . . . . . . . . . . 47 3.8 Numerical solution to Sod’s problem (3.18). The total number of elements K and polynomial order N is varied. O’s: Q-ENO recon- struction. X’s: DK reconstruction. Solid line: exact solution. . . . . . 48 3.9 Numerical solution to the Lax problem (3.19). K = 200 elements and N = 2 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. . . . . . . . . . . . . . . . . 50 3.10 Numerical solution to the Lax problem (3.19). A fixed value of K = 200 grid points with varying N = 1, 2, 3, 4 order is given. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. . . 51 3.11 Numerical solution to the Lax problem (3.19). Solutions of order N = 1 and varying K = 100, 200, 400, 800 are given. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. . . 52 3.12 Numerical solution to the Shock density wave problem (3.20). K = 400 elements and N = 3 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted (one point per element). O’s: Q- ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 54 3.13 Numerical solution to the Shock density wave problem (3.20). K = 400 elements with varying order N = 2, 3, 4, 5 order are used. One point per element is plotted. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. . . . . . . . . . . . . . . . . 55 xiii 3.14 Sequence of refined grids with K = 146, 568, 2310 (unstructured) tri- angular elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.1 Illustration of the formation and coalescence of melt channels. Adi- abatic decompressional melting occurs within the mantle. Channels migrate upwards, with impenetrable solid barriers guiding the melt towards the mid-ocean ridge. Picture taken from [31]. . . . . . . . . . 60 4.2 The Mid-Atlantic Ridge, a 10,000 km long string of volcanoes. Colors indicate relative age of the crust under the ocean. Red is youngest, with older ages as one approaches the continents. Picture taken from [31]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.3 Base states p¯, c¯f for the difference choice of boundary conditions. Top: p(0) = cf (0) = 0. Bottom: modified boundary conditions as given in Eq. (4.58). The domain for cf is shortened to z ∈ [0, 0.02], as outside of this it is essentially parallel to cf = ceq f = z. . . . . . . . . . . . . . 74 4.4 Gauss-Lobatto quadrature nodes for order N = 20. The nodes are given over the one-dimensional reference element ξ ∈ [−1, 1]. . . . . . 79 5.1 Linear stability analysis, as computed in Eqs. (5.2)-(5.10). γ = 0.00025. Top Left: computed growth rates σ vs. perturbation fre- quency k, for varying Da. Top Right: convergence in φ˜f , p˜, σ as Da → ∞. A fixed perturbation frequency of k = 100 is chosen. Bottom left/right: φ˜f /˜ p eigenfunctions for k ∈ [0, 240], plotted in increments of 20. All eigenfunctions have been normalized so that φ˜f (zH ) = p˜(zH ) = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 xiv 5.2 Comparison of predicted growth rates σ from linear stability analysis (LSA) from Eqs. (5.2)-(5.10) to those actually measured from the transient simulations (4.38)-(4.43). Three perturbation frequencies, k = 5, 150, 225 are examined for the transient simulations. Growth rates in the transient simulations are computed by considering the rate of change of the magnitude of the Fourier coefficient correspond- ing to frequency k. Top: Equilibrium formulation. Bottom: Kinetic formulation (Da = 106 ) with surface area coefficient l = 0.5. For both cases we take m = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.1 Growth rates σ vs. horizontal perturbation frequency k from linear stability analysis. γ = 0.00025. Left: fixing l = 0, we vary the porosity component φm f in bulk viscosity. A reaction rate constant of Da = 106 is used, which is essentially converged to the equilibrium case. Right: fixing m = 0, we vary the surface area component φlopx in reaction rate. A reaction rate constant of Da = 104 is used. Taking both l, m → 1 stabilizes the system. . . . . . . . . . . . . . . . . . . . 97 6.2 Comparison of input (bottom) perturbation to output (top) channels formed. Porosity φf is shown for a reaction rate of Da = 10 up to a dimensionless time of t = 20. Porosity perturbations from Eq. (6.17) are used: For the first 7.33 time units, a low-frequency perturbation with noise in the first five Fourier modes is used. For 7.33 < t < 11, a higher frequency perturbation, with noise in the first 10 Fourier modes, is given. For 11 < t < 20, a regulated series of three Gaussian bumps is prescribed, with the middle disappearing every other time unit. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 xv 6.3 Transient simulation with perturbed φf inflow as pictured in the bot- tom panel of Figure 6.2. For this case, however, we take Da = 1000. The solution is examined at time t = 2.8, which corresponds to the perturbation in the first time regime in Eq. (6.17). Five channels of melt are observed, corresponding to the perturbation noise at in- flow contained with the first five Fourier modes. Melt channels in φf correspond one-to-one with dissolved regions of φopx . . . . . . . . . . 104 6.4 Solution φf and φopx at outflow (z = 3.66) for varying Da. The same time t = 2.8 as in Figure 6.3 is chosen, with perturbation defined in as in Eq. (6.17) and plotted in Figure 6.2. Higher Da produces deeper channels of melt and opx dissolution. . . . . . . . . . . . . . . . . . . 106 6.5 Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Top row: Two-dimensional profiles for Da = 500. Middle and bottom rows: solutions at in the channel mid- point x = 0.5 as a function of z. Once the opx mineral is exhausted, the porosity drops and compaction rate rises sharply. . . . . . . . . . 108 6.6 Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Reaction constant is Da = 500. Soon after the opx mineral is exhausted, the melt channel bifurcates into a “Y” formation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.7 Top: fluid velocity vector Vf = (vx , vz ), for profiles presented in Fig- ure 6.6. For vx , negative (positive) values represent leftwards (right- wards) flow. Before the channel bifurcates, the fluid tends to flow more inwards towards the channel, while after the bifurcation it flows outwards. The vertical component vz resembles the porosity. Bot- tom: fluid velocity streamlines. The fluid tends to avoid low-porosity compaction regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 xvi 6.8 Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Reaction constant is Da = ∞. Nu- merical results are very similar to Figure 6.6, but are available at considerably less computational expense. . . . . . . . . . . . . . . . . 112 6.9 Comparison of varying Da behaviors at solution outflow z ≈ 3.66, time t = 2, for transient solution with fixed boundary perturbation given in Eq. (6.18). The equilibrium solution closely resembles the kinetic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 xvii Chapter 1 Discontinuous Galerkin Methods The two projects presented in this thesis build upon the foundation of the discon- tinuous Galerkin (DG) finite element method. The method will be presented in this chapter, with discussion given to implementation-oriented concerns such as modal and nodal representations and stability for nonlinear problems. A brief discussion will be given for explicit Runge-Kutta time integration methods. Finally, we will discuss the lowest-order implementation of DG, the finite volume method. 1.1 History The discontinuous Galerkin method is a popularly emerging choice for the numerical solution to partial differential equations. The first DG method was introduced in 1973 [47] as a means of solving the steady-state neutron transport equation. Major developments in the DG method were introduced in [12] with the combination of explicit, nonlinearly stable Runge-Kutta methods [50] with an improved slope limiter [49], allowing formal accuracy in regions where the solution is smooth, sharp shock resolution, and convergence to the entropy solution. Further progress was made with the extension to high-order accuracy in [13] and one-dimensional systems in [11]. In [10], these results were extended for multidimensional scalar problems, and finally 1 2 to systems in [14] with the introduction of a special class of slope limiters. 1.2 Motivation DG combines many attractive features of the classical finite element (FE), finite vol- ume (FV), and finite difference (FD) methods. In short, DG schemes, with a poly- nomial solution and test function basis, appear similar in presentation to continuous finite element method (“Galerkin”). However, there are no continuity requirements on the solution at element interfaces (“discontinuous”), which more resembles the nature of the finite volume method. Indeed, the lowest order DG scheme is equivalent to FV, as the evolved degree of freedom is a cell average. Generally, however, one chooses to implement DG for its ability to provide a high- order solution, which contrasts with the inherently low-order FV, and its freedom to be implemented on unstructured grids, which contrasts with FD. DG is ideally suited for conservation laws and problems with propagating waves, as the discontinuous solution representation at cell interfaces decouples the global solution into a collection of local ones. Therefore, rather than having large global mass and stiffness matrix operators, we have local operators that may be tailored specifically to the physics of the problem at hand. This is in contrast to the continuous FE method. One obstacle that DG faces, however, is the well-known order barrier theorem from Godunov which states that high-order schemes are necessarily non-monotonic. For smooth solutions this poses no problems, but in the presence of strong discon- tinuities, the numerical solution will contain non-physical oscillations. Techniques such as slope limiting have been used to resolve this, but generally require fine- tuning of solution and/or mesh-dependent parameters in order to avoid reducing the order of the scheme in smooth regions. More robust limiter choices such as WENO [46, 44] and Hermite WENO [43, 45] shows great promise, but are less suitable for unstructured grids. In Chapter 3 a hybrid scheme, which switches between the dis- 3 continuous Galerkin method and a finite volume method with high-order function reconstruction for the evaluation of numerical fluxes, is presented. In this work posi- tive research results (high-order, non-oscillatory) are given for the numerical solution to one-dimensional conservation laws on unstructured grids. 1.3 Basic Formulation In this section we present the basic formulation of the discontinuous Galerkin method scheme for scalar conservation laws. These are generally written as ∂u (x, t) + ∇ · f (u(x, t), x, t) = 0, x ∈ Ω ⊂ Rn ∂t u (x, t) = g (x, t) , x ∈ ∂Ω u (x, 0) = u0 (x), x ∈ Ω (1.1) where the flux function f is given, and boundary conditions g are available to ensure well-posedness of the PDE. T k , locally, the N -order DG solu- S Tesselating the domain into elements Ω = k tion to (1.1) is chosen from the finite polynomial basis Vhk = φk : φk ∈ P N (T k ), k = 1, 2, ..., K ,  (1.2) where P N (T k ) is the set of polynomials of degree N with support in element T k . We define the dimension of Vhk to be Np , and for general tetrahedral elements in n-dimensional space, the scaling between N and Np is   n+N n Np =  ∼ N (1.3) N n! Returning to (1.2), the basis Vhk is selected as both the solution basis and the 4 test function basis, as is standard for a Galerkin scheme. There are no continuity requirements on the solution across element boundaries ∂T k . For simplicity we as- sume that the same basis Vhk is used for each element, although we may certainly vary this, i.e., adaptively choosing N based upon the expected smoothness of the so- lution. Locally, then, the solution may be represented as either a linear combination of basis functions Np X ukh (x, t) = uˆkj (t)φkj (x), x ∈ T k (1.4) j=1 or equivalently as a linear combination of nodal values Np X ukh (x, t) = ukh (xkj , t)ljk (x), x ∈ T k . (1.5) j=1  Np Here, a set of gridpoints xkj j=1 has been introduced along with the corresponding Lagrange interpolating polynomials ljk . For general tetrahedral mesh, any two elements may be related to one another via a linear transformation, and so for implementation purposes we may perform all computations on a reference element T E with interpolation points ξj and then map back to T k . Thus, all operators are local, with only the need to store the Jacobians of each element’s transformation to T E . On the reference element, then, the two solution representations (1.4) and (1.5) become (suppressing time dependence) Np X uE h (ξ) = uˆj φj (ξ) j=1 Np X = u (ξj ) lj (ξ) (1.6) j=1 5 The degrees of freedom uˆj and u (ξj ) are linked via the Vandermonde matrix Vu ˆ=u (1.7) where Vij = φj−1 (ξi ) , u ˆ i = uˆi , ui = u (ξi ) . (1.8) Practically speaking, it is easier to perform mathematics on the modal repre- sentation (1.4) if one chooses orthonormal φkj . However, it is easier to visualize the scheme (and thus, the solution) by considering the nodal representation (1.5). The Vandermonde matrix is typically found embedded in other operators in the scheme (e.g. the mass and stiffness matrices), so that the solution remains in nodal form while computations are performed in modal form. The choice of which basis functions φkj or which interpolation values ξjk to use are of significance when one considers issues such as matrix conditioning; see [26, 59] and the references therein for more details. Typically for φkj we use the orthogonal Legendre polynomials, a special subset of the Jacobi polynomials, with corresponding Gauss-Lobatto interpolation points ξjk (See Sec. 4.4 for more discussion). In the latter case the basis functions become the Lagrange interpolating polynomials ljk . Greater discussion is given to the derivation of the Gauss-Lobatto points in Sec. 4.4.2. Having defined the solution space, we proceed with defining the scheme by consid- ering the variational form of (1.1). Multiplication by test function ljk and integration over T k yields   R ∂ukh Tk ∂t + ∇ · fhk ljk dx = 0, j = 1, 2, ..., Np (1.9) where the discretizations ukh , fhk have now been introduced as approximations to u and f on element k. Using the divergence theorem we then integrate by parts to 6 arrive at the residual statement   R ∂ukh k H ∗ Tk l ∂t j ˆ · fhk ljk dx, j = 1, 2, ..., Np . − fhk · ∇ljk dx = − ∂T k n (1.10) Alternatively, we may integrate by parts again and have   R ∂ukh H  ∗  Tk ∂t ˆ · fhk − fhk ljk dx, j = 1, 2, ..., Np . − ∇ · fhk ljk dx = − ∂T k n (1.11) Equations (1.10) and (1.11) are the weak and strong forms, respectively, of the discontinuous Galerkin method. The two forms are mathematically equivalent but differ in implementation. Noting that the numerical solution ukh is discontinuous ∗ across element boundaries ∂T k , we must introduce a new approximation fhk , the numerical flux, to compute the right hand side of these equations. Whereas the fhk terms and its derivatives are computed exclusively on the element T k , the numerical ∗ flux fhk depends on values taken both from T k and its neighbors. As a result, the local solutions become coupled, allowing for the global transfer of information. Further discussion concerning the numerical flux can be found in Section 1.6.1. In Sec.4.4.4 an implementation of the strong form is given for a system of conservation laws. Having now arrived at a spatial discretization (1.10) or (1.11), we may integrate the solution in time with a suitable time-stepping scheme. Next, we discuss the explicit Runge-Kutta method, a special class of such schemes. 1.4 Explicit Runge-Kutta Time Integration The equations resulting from (1.10) and (1.11) reflect a method of lines approach, where the spatial discretization of (1.1) reduces the partial differential equation to a system of ordinary differential equations. This results in a semi-discrete scheme of 7 the form duh = Lh (uh , t) (1.12) dt for some corresponding operator Lh . We then advance the solution with a suitable time integration scheme. One such method comes from the class of explicit Runge- Kutta (ERK) schemes. General m−order ERK schemes [52] advance the solution unh = uh (·, tn ) in time to un+1 h = uh (·, tn+1 ) = uh (·, tn + ∆t) by evaluating the operator Lh in m stages, v (0) = unh i−1 X v (i) = αi,k v (k) + ∆tβi,k Lh v (k) , tn + ∆t(k)  , i = 1, ..., m, k=0 v (m) = un+1 . (1.13) The m = 2, 3 and 4 versions of the classic ERK (1.13) are popular choices for time integration, but come with the drawback of needing to store each of the intermediate stages v (i) . An alternative to the fourth order scheme is the low-storage version LSERK [7], p(0) = unh ,   v (i) = a v (i−1) + ∆tL p(i−1) , tn + c ∆t , i h i i ∈ [1, ..., 5] : (1.14)  p(i) = p(i−1) + b v (i) i un+1 h = p(5) Although the fourth-order LSERK requires an additional stage of computation, there is only one additional storage level v (i) , thereby significantly decreasing the storage cost. Moreover, LSERK potentially allows for a larger stable timestep ∆t than its classic ERK counterpart; see [28] for details. 8 Runge-Kutta methods are the most popular selection for time-integration. For more choices, see [6, 23, 24]. 1.5 Problems With Nonlinear Problems Having now walked through a basic discontinuous Galerkin scheme, we turn to focus on one potential problem that may arise with a nodal implementation. It can be proven [29] that the discontinuous Galerkin solution to (1.1) converges to the entropy solution if f is convex. Practically speaking, however, caution should be urged in implementation. If the solution is represented by storing its nodal values as in (1.5), then we can evaluate the flux as Np X fhk = f (ukh (xkj , t))ljk (x), (1.15) j=1 that is, by simply reading uh at the grid points and evaluating for fh . This is an easy and straightforward way to compute the flux, but only accurate if f itself is able to represented by the basis chosen to represent uh . For example, if f (u) = u4 , uh ∈ P N , then fh ∈ P 4N , and the representation (1.15) will have insufficient number of gridpoints to completely represent the solution. Thus, we would incur aliasing errors that could potentially lead to instability in the numerical solution R when evaluating stiffness terms such as the integral T k (∇ · fh ) ljk dx. Obviously, if we were to increase the order of approximation N , the increased resolution would eventually be able to sufficiently capture fh . While this is a robust way to solve the problem it will certainly increase the computational cost and often is not practical. Further discussion of this can be found in [22, 27]. A different solution is to add artificial dissipation to the problem as a means of filtering out the unresolvable, highest-order modes of the solution. Typically, this involves solving a modified form of (1.1) in which we add additional terms of the form 9 ∂ 2s u ∂x2s for some (even) integer 2s. Instead of worrying about the operators involved with these terms and how they would be incorporated into the scheme, we may simply post-process the solution to achieve the goal at hand: Np   X j−1 uh (ξ) → F uh (ξ) = σ uˆj φj (1.16) j=1 N where the exponential filtering function σ is defined as  Nc  1, 0 ≤ η ≤ ηc = N σ(η) = (1.17)  exp (−α ((η − η ) / (1 − η ))s ) , η ≤ η ≤ 1. c c c From (1.16)-(1.17) we can see that the highest order modes in the solution are damp- ened, depending upon the parameters chosen. The parameter Nc acts as a cutoff frequency whereby lower-order modes are unfiltered. α is a dissipation parameter, usually taken to be on the order of α = − log (M ), where M is machine precision. Whereas α tends to remain fixed, the strength parameter s adds a dimension of vari- ability; a good starting point for this value is s = 16. Higher Nc and s correspond to less filtering. For the second project in this thesis presented in Chapter 4, which presents high- order methods for the reaction infiltration instability, filtering is used to dissipate the aliasing errors introduced from the nonlinearity in the equations. The numerical schemes in this work implement a Fourier spectral method; the exponential filter presented may also be applied in this context as well. Details of the implementation can be found in [27]. 1.6 Introduction to Finite Volume Schemes In the lowest order (N = 0) setting, the discontinuous Galerkin method reduces to a finite volume scheme. In this section we will step through a one-dimensional imple- 10 mentation, providing further discussion of the numerical flux which was previously introduced in the general formulation of the discontinuous Galerkin method. We consider the one-dimensional scalar conservation law ∂u(x, t) ∂f (u) + = 0, x ∈ Ω ⊂ R, t > 0 ∂t ∂x u(x, 0) = u0 (x) (1.18) subject to suitable boundary conditions so as to guarantee well-posedness of the problem. We assume a domain Ω = K k k S k=1 T , where the elements T are the intervals T k = xk−1/2 , xk+1/2 ,   (1.19) so that T k ≡ xk+1/2 − xk−1/2 . Using this notation, we may denote f on cell  boundaries, f k±1/2 = f xk±1/2 . The finite volume scheme is derived by integrating (1.18) over each interval T k , uk d¯ 1 = − k f k+1/2 − f k−1/2  (1.20) dt |T | where we have introduced the notation for the cell average u¯k , Z xk+1/2 k 1 u¯ (t) ≡ k u(z, t)dz. (1.21) |T | xk−1/2 These averages are the fundamental degrees of freedom to be solved for. The solution is therefore discontinuous across cell interfaces xk±1/2 , and so we must approximate the flux f in these locations. This approximation, the so-called numerical flux is denoted ∗   f k+1/2 = h (u−)k+1/2 , (u+)k+1/2 (1.22) where the notation (u−)k+1/2 and (u+)k+1/2 refers to left- and right-sided approxi- mations of u, with respect to the interface between elements T k and T k+1 . 11 1.6.1 Numerical Flux The numerical flux function h = h(a, b) introduced in Eq. (1.22) must satisfy the following three properties: • h is Lipschitz continuous with respect to both arguments • h is nondecreasing with respect to its first argument and nonincreasing with respect to its second argument • h is consistent with the original flux f , that is, h(a, a) = f (a) There are a great number of numerical flux functions that may suffice for a given problem. A comprehensive listing of these can be found in [36]. Perhaps the easiest and most stable (dissipative) choice is the Lax-Friedrichs flux 1 hLF (a, b) = [f (a) + f (b) − C (b − a)] (1.23) 2 where the constant C is df C = max . (1.24) u du The maximum may be applied over the entire domain Ω (global) or in a more re- stricted range (local), such as in small neighborhood of the element interface. In this case C will vary throughout the domain. The global Lax Friedrichs flux provides a more dissipative scheme. An alternative approach to (1.23) adds a control 0 ≤ α ≤ 1 to the dissipative term by taking 1 h(a, b) = [f (a) + f (b) − C (1 − α) (b − a)] . (1.25) 2 In this way, choosing α = 0 recovers (1.23) while choosing α = 1 results in a purely central flux with no dissipation. The latter is the preferred choice when dealing with problems with no preferred direction of propagation, e.g. the diffusion equation. 12 Returning to the numerical scheme (1.20), then, we require a reconstruction of ∗ the values u(xk±1/2 ) to approximate f k±1/2 . Once this is known we may evolve the solution in time utilizing an ERK method, as presented in Section 1.4. 1.7 Conclusion This chapter has introduced the discontinuous Galerkin finite element scheme, with additional discussion of its lowest-order (N = 0) version, the finite volume scheme. Essential for both formulations is the numerical flux. For the N > 0 discontinuous Galerkin method the numerical flux may be readily computed by considering the solution values on all sides of the interface in question. For the finite volume scheme, however, this information is not available. In the next chapter we will detail solution reconstruction algorithms which provide the needed quantities. Part I Hybrid Discontinuous Galerkin Methods 13 Chapter 2 ENO-Type Reconstructions In both the discontinuous Galerkin method as well as the finite volume method, the discontinuous solution is connected to its neighbors via the numerical flux. The nu- merical flux is computed via an approximate Riemann solver, considering the values of the multiply-defined solution across element boundaries. For the discontinuous Galerkin method of order 1 and greater, these values are naturally available due to the polynomial basis which represents a continuous local expression of the solution. However, for finite volume schemes the fundamental degrees of freedom used are merely the cell-averages. The computation of the numerical flux, therefore, requires a reconstruction procedure whereby one uses these degrees of freedom to compute an approximation of the solution at element boundaries. ENO, or essentially non-oscillatory methods, first presented in [25], may be used to provide an oscillation-free reconstruction capable of preserving high-order features. The main idea of ENO is to pick the “best” solution reconstruction among those computed from a collection differing stencils. WENO [37, 30] methods improve upon this idea by instead considering a weighted sum of reconstructions across all stencils, thereby improving the overall order of accuracy. The particular weights to use are computed via a nonlinear smoothness indicating functional of the cell- averaged degrees of freedom, with nonsmooth reconstructions receiving weights on 14 15 the order of the truncation error in order to attain the ENO property. The main difficulty involved in the reconstruction process is that of maintaining accuracy. In particular, high-order features of a solution may be incorrectly com- puted without a suitably high-order reconstruction. At the other end of the spec- trum, solutions with strong shocks may suffer non-physical oscillations or excessive smearing due to poor reconstruction. Moreover, standard WENO schemes [37, 51] are inherently one-dimensional in nature, performing reconstructions direction-by- direction on structured mesh. For the more complicated case of an unstructured mesh in multiple dimensions, ENO-type reconstructions are hindered by the need for potentially large stencils and complicated stencil-searching algorithms to detect the best stencil(s) from which to build the solution. Recent developments have in- cluded Hermite WENO reconstructions [43, 45], which include both the cell average and first derivatives as the degrees of freedom. These show promise in reducing the stencil size, but still carry the requirement of a structured mesh. A Hermite WENO scheme for unstructured mesh is given in [38], but is as yet only presented for linear reconstructions. In this chapter, we will first discuss the general attributes of ENO-type recon- structions. Next, we will detail the implementation of two particular high-order ENO-type reconstruction algorithms, QENO and DK, proposed in related forms in [41] and [17, 18], respectively, will be discussed. The former resembles an ENO scheme by reconstructing a single polynomial while the latter is a WENO scheme, with multiple reconstructions weighted in a nonlinear fashion. Both involve solving least-squares problems, which allows for grace in the stencil-searching procedure to seek out smooth data from which to reconstruct the solution. 16 2.1 Properties of the ENO Class The function reconstruction procedure is as follows: given a stencil of cell averages S¯k = u¯k−p , u¯k−p+1 , ..., u¯k+q−1 , u¯k+q ,  (2.1) we wish to build an approximating polynomial used to evaluate u on cell boundaries, uk±1/2 . Perhaps the simplest way to do this would be via an averaging, i.e. 1 k uk+1/2 = u¯ + u¯k+1/2 .  (2.2) 2 With larger stencils, one could expand this interpolation to higher order accuracy. Unfortunately, as is well-known in function interpolation theory, such reconstruc- tions of order two or greater will necessarily produce oscillatory data should the stencil include discontinuous solutions. ENO-type reconstructions are a special class of reconstruction algorithms that produce high-order, (essentially) non-oscillatory reconstructions. We seek a reconstruction of this class with the following properties: • Conservation of mean: We must preserve the fundamental degrees of free- dom (1.21). • Accurate representation of smooth functions: To maintain accuracy for smooth functions it is required that for any smooth function u(x), the corresponding local reconstruction Ri (x − xi ) will satisfy Ri (x − xi ) = u(x) + O(∆xN +1 ). (2.3) • Compact support: One very attractive property of the discontinuous Galerkin scheme is that locally, the solution only depends upon its immediate neighbors. Unfortunately this is not necessarily the case with finite volume schemes as the reconstruction of function values at cell boundaries may involve the need to 17 consider values from a wide stencil of neighbors in addition to the local element itself. The construction of a one-dimensional polynomial of (high-order) degree N , for example, requires at least N + 1 degrees of freedom. How these degrees of freedom are chosen will determine the stencil. • ENO property: Godunov’s theorem guarantees that linear schemes of order greater than one necessarily will produce non-physical oscillations. One may circumvent this result, then, by introducing nonlinearity into the scheme. That is, the reconstruction itself depends upon a nonlinear functional of the stencil averages (2.1). 2.1.1 Smoothness Indicators The essentially non-oscillatory, or ENO, property is introduced to limit non-physical numerical oscillations in regions with steep gradients. ENO reconstruction schemes assume that a smooth polynomial exists locally even in regions with shocks, if only we use the appropriate selection (stencil) of data to construct it. Weighted ENO, or WENO, advances this by applying nonlinear weights to combine all reconstructions into a higher-order, essentially non-oscillatory polynomial. Central to the computa- tion of these weights is the smoothness indicating functional N −1 Z xk+1/2 2 dl pr  X k 2l−1 SI r = T dx, (2.4) l=1 xk−1/2 dxl where pr is the rth interpolated polynomial of degree N − 1. SI r is then just a sum of squares of L2 norms of pr (i.e., a Sobolev norm) over the cell T k . The T k scaling removes mesh dependency in the functional. Assuming a uniform mesh, SI r reduces to a simple nonlinear function of cell averages. WENO proves to be quite robust, but is unfortunately restricted to uniform mesh. We will present two different ENO-type reconstructions, [41] and [17, 18], that do not carry such a restriction. The main distinguishing factors among these and other ENO-type schemes are the smoothness 18 indicating functions and the weights used in combining the reconstructions. 2.2 QENO Reconstruction Scheme In [41], an algorithm is proposed which uses data-dependent non-linear weights to construct a single polynomial. The self-described “Quasi-ENO” scheme uses a least- squares approach to reconstruct a high-order, single polynomial which preserves the average over the native element. Nonlinear, data-dependent weights are used to assign priority to nearest-neighbors’ information (cell average) in the least-squares solve. A slightly modified version of this reconstruction is presented below, and is implemented in the Hybrid scheme. We will refer to this scheme as QENO. 2.2.1 Conservation of Mean To preserve the mean, the reconstruction Ri of an N th -order polynomial on cell i is N  1 ∂ k u  X k k  Ri (x − xi ) = u¯i + (x − xi ) − x¯i (2.5) k=1 k! ∂xk i where the geometric coefficients are given as Z 1 x¯ki = (x − xi )k dx (2.6) |Vi | Vi By inspection it is seen that Eq. (2.5) preserves the mean. 2.2.2 Stencil Selection For each cell i there is an associated stencil S i = {i − p, i − p + 1, ..., i + q − 1, i + q} 19 where p ≥ 0 and q ≥ 0 are the numbers of neighbors to include in the reconstruction on the left and and right, respectively. For implementation purposes, as will be seen later, the native cell i is not included in the stencil index S i , as its ”information” (namely, the cell average u¯i ) is used to enforce conservation of mean and will not be included as an equation in the least-squares problem. For this research the values of p = q = N are taken, thereby giving a total stencil of size 2N . That is, there will be double the amount of information given than what is required to reconstruct a degree-N polynomial. Assuming a simple ordering of grid points (xi−1 < xi < xi+1 for all i), a volume from stencil S i and its corresponding function average are denoted Sji and u¯ij , respectively, where 1 ≤ j ≤ 2N and   i + j − N − 1, 1≤j≤N  i Sj = ij = . (2.7)  i+j−N N + 1 ≤ j ≤ 2N  Appropriate modifications are needed near boundaries. 2.2.3 Reconstruction Algorithm To accomplish exact reconstruction of degree N polynomials, a least-squares proce- n k oN dure is used to compute the N derivative terms k!1 ∂∂xuk in the expansion of Ri i k=1 in Eq. (2.5). These derivatives will be computed by minimizing the predicted mean  2N value of Ri over all control volumes u¯ij j=1 in the stencil S i . Each of these mean values can be expressed explicitly as    N  1 ∂ k u  1 Z  Z  1 X k k  Ri (x − xi )dx = u¯i + (x − x ) dx − x ¯ (2.8)  i i Vi  k! ∂xk i Vij  j  Vij k=1  Vij  20 To simplify this expression, we extract the geometric terms and define Z 1 xˆkij = (x − xi )k dx − x¯ki . (2.9) V ij Vij Thus, (2.8) is rewritten as N 1 ∂ k u k Z 1 X Ri (x − xi )dx = u¯i + xˆ (2.10) k! ∂xk i ij V i j k=1 Vij For neighbor Vij in the stencil of Vi , the difference in the predicted mean value and the actual value is N  1 ∂ k u k Z  1 X  Ri (x − xi )dx − u¯ij = x ˆ ij + u ¯ i − u ¯ i j (2.11) V i j k=1 k! ∂xk i Vij The least-squares problem is tentatively formed, then, by setting Eq. (2.11) equal to zero for all volumes in the stencil. Solving for the derivatives, this is a system of N unknowns with 2N equations coming from each of the volumes in the stencil. The least-squares solution to such a system, however, enforces Eq. (2.10) equally for all neighbors in the stencil. This is clearly unacceptable, as one would like to bias the influence of neighbors in such a way that 1. preference is given to volumes physically closer to the reconstruction volume (compact support), and 2. the contribution of discontinuous data in the reconstruction is reduced to the level of truncation error (ENO property), yet 3. the contribution of smoothly connected data is unperturbed (high-order accu- racy) For reconstruction on volume Vi these items can be accomplished by appropriately 21 selecting weights wij that will weight each of the j = 1, ..., 2N equations from (2.10) in the least squares problem. 2.2.4 Weights and the ENO Property The weights will have two components: geometric, which can be computed during a preprocessing stage, and a smoothness indicator, which is nonlinearly data-dependent and must be computed for each reconstruction. The geometric contribution is simply g 1 wij = (2.12) |xj − xi |2 The data-dependent component is more complex. It is presented as , ! u¯j − u¯i N +1 DD wij =1 ε + (2.13) xj − xi where ε = 10−10 is chosen to prevent division by zero. The divided difference u¯j −¯ui N +1 xj −xi is the smoothness indicator. Previous work [41] solved the least squares solution twice: once with only the geometric component of the weights (i.e. a data-independent problem), and then again with the same weights as above, except with the residual of the first results multiplied into the divided difference in the denominator of wDD . However, with the assumption that the Hybrid will only invoke the L2 -ENO scheme in regions with questionable smoothness, it is unnecessary to use the residual component as by circumstance it will automatically be significantly larger than ε and therefore scale each row of the least-squares problem equally. The weights then can be formally defined as !−1 u¯i − u¯ij N +1 g DD 1 wij = wij wij = · ε + (2.14) xi − xi 2 xi − xi j j 22 Thus, the least-squares problem on control volume i with neighbors indexed j = 1, ..., 2N is the system Ld = b, (2.15) where Ljk = wij xˆkij 1 ∂j u dj = j! ∂xj (2.16) i bj uij − u¯i ). = wij (¯ 2.3 DK Reconstruction Scheme The reconstruction algorithm presented in Sec. 2.2 constructs a single polynomial by solving a least-squares problem with different weights corresponding to the relative importance of preserving neighboring cell averages. An alternative idea follows the more traditional WENO scheme: compute multiple polynomial reconstructions over differing stencils, and then combine these together in a non-linear fashion to arrive at a final high-order, non-oscillatory reconstruction. The scheme presented is related to the procedures of [18] and [17], and is named DK after the authors of this work. 2.3.1 Polynomial Representation One can view the fundamental degrees of freedom to be reconstructed as either nodal values or modal coefficients of the polynomial. These forms are N X N X u(x) = u(xj )lj (x) = wˆj Pj (x), (2.17) j=0 j=0 respectively. The DK scheme reconstructs the latter. 23 2.3.2 Stencil Selection The domain is tessellated into conforming elements T (m) , K [ Ω= T (m) . m=1 On each of element T (m) a selection of ns reconstruction stencils will be selected by considering collections of ne neighboring elements, ne [ (m) Sj = T j(k) , 1 ≤ j ≤ ns (2.18) k=1 By convention we take j(k) = m for all stencils, that is, the first element in each stencil is the local element itself. This is necessary for the conservation of the cell average. More precisely, the mapping function j(k) carries a dependence on the local element m, and so can be thought of as j = j m (k). We will consider now an arbitrary element and so suppress this dependence on m in the notation. The problem that remains is to select fj (k). For the one-dimensional case the stencil selection is straightforward. To recon- struct a polynomial or degree N we need at least N + 1 cell averages, and therefore we must have ne ≥ N + 1. Each stencil will have the form (m) Sj = {m, m − r, m − r + 1, ..., m − 1, m + 1, ..., m + s − 1, m + s} (2.19) By convention the first stencil (j = 1) is a central stencil, chosen to contain the near- est neighbors. For this case one may take r = s = b (N + 1)/2c, the greatest integer less than (N + 1)/2. Note that this will in some cases lead to an overdetermined system. This is acceptable, and even necessary in higher dimensions. For the second and third stencils (j = 2, 3) we may pick a left and right bias, 24 respectively, leading to r = N, s = −1 and r = −1, s = N . Depending on whether N is odd or even, and thus whether the the first stencil is of length N + 1 or N + 2, it may be computationally convenient to also expand the second and third stencils to include an extra element. Thus one might instead choose r = N + 1, s = −1 and r = −1, s = N + 1. Stencil selections for higher dimensions are a bit more complicated. For two- dimensional triangular elements, there are ns = 7 stencils: a central stencil, three primary stencils in each of the directions normal to the sides of the element, and three additional “reverse WENO” stencils. The latter three are computed by con- sidering instead the vertices of the triangle and then taking the elements spanned by the reverse (negative) directions of the correspodning primary stencils. In three dimensions, the number of stencils is ns = 9. For a more thorough discussion of the stencil-building process, see [17]. 2.3.3 Reconstruction Algorithm Polynomial reconstructions R(m) (x) on element T (m) will be computed on each stencil by requiring conservation over each cell in the stencil. That is, Z 1 (m) R(m) (x)dx = u¯j(k) , 1 ≤ k ≤ ne , j(k) ∈ Sj (2.20) |T (j(k)) | T (j(k)) It is helpful for the sake of implementation to address how this integral will be computed. For numerical quadrature we would like to transform coordinates so that all integrals on a reference element T E . This transformation is denoted as x = x(T (m) , ξ), ξ = ξ(T (m) , x) (2.21) The Jacobian of this transformation is ∂xi Jij = (2.22) ∂ξj 25 and its determinant is |J|. Therefore (2.20) may be written as |J| Z (m) R(m) (ξ)dξ = u¯j(k) , 1 ≤ k ≤ ne , j(k) ∈ Sj (2.23) |T (j(k)) | T˜(j(k)) The integrals in Eq. (2.23) are over the (j(k))-indexed elements in the ξ-coordinate space; still, for the sake of implementation, these will be each transformed again, ˜ this time to the reference element T E in ξ-coordinate space. This transformation is denoted as ˜ ξ˜ = ξ( ξ = ξ(T˜(j(k)) , ξ), ˜ T˜(j(k)) , ξ) (2.24) ˜ Now, Similar to Eq. (2.22), the Jacobian of this second transformation is denoted J. the element to be integrated over is simply the reference element T E and so the conservation equations (2.23) can be written as |J||J|˜ Z ˜ ξ˜ = u¯j(k) , 1 ≤ k ≤ ne , j(k) ∈ S (m) R(m) (ξ(T˜(j(k)) , ξ))d (2.25) j |T (j(k)) | TE Assuming the case of straight-edged elements (e.g. triangles and tetrahedrons), the Jacobian J is constant and we have the relationship |T (j(k)) | |J| = . (2.26) |T˜(j(k)) | Likewise, we have ˜ = |T˜(j(k)) | |J| (2.27) |T E | Thus, the Jacobian factors cancel out and, using the convention of summation over repeated indices and recalling the representation from Eq. (2.17), the conservations equations may be written  Z  1 ˜(j(k)) ˜ ξ˜ wˆl = u¯j(k) , 1 ≤ k ≤ ne , j(k) ∈ S (m) Ψl (ξ(T , ξ))d j (2.28) |T E | TE 26 Using the convention of summation over repeated indices, Eq. (2.28) can be suc- cinctly written in tensor notation as Akl wˆl = u¯k (2.29) where  Z  1 ˜ (j(k)) ˜ ˜ Akl = Ψl (ξ(T , ξ))dξ wˆl (2.30) |T E | TE and u¯k = u¯j(k) (2.31) We finally have a linear system to solve which will produce a reconstruction for stencil j, 1 ≤ j ≤ ns . If we denote the number of degrees of freedom to be solved for as Np (N + 1 in one dimension and 21 (N + 1)(N + 2) in two dimensions), the number of neighbors ne (including the native element itself) needed for reconstruction must be greater than or equal to Np . In one dimension we may take ne = Np or ne = Np + 1, depending on whether one would like a balanced, central stencil (that is, equal number of left- and right- neighbors). However for higher dimensions a greater stencil size should be included, so as to handle any potential singularities in the linear system that may arise due to the unstructured mesh. In [17] the values ne = 1.5Np and ne = 2Np are suggested for two and three dimensions, respectively. Ensuring Conservation of Mean The linear system (2.29) will become overdetermined when ne > Np , and so an additional constraint must be added to ensure each reconstruction is mean-preserving over its native element T (m) . Recall that the function representation is Ψl wˆl , where Ψl are the basis functions and wˆl are the degrees of freedom. If we use Legendre polynomials as basis functions, the necessary constraint is 27 1 j(1) 1 wˆ0 = u¯ = u¯(m) (2.32) γ0 γ0 where γ0 is the normalization constant, dependent upon whether one is using nor- √ malized (γ0 = 2) or non-normalized (γ0 = 1) basis functions. This constraint is written in tensor notation as 1 Cl wˆl = Ri u¯i (2.33) γ0 1 where all entries of Cl and Ri are zero except for their first, C0 = γ0 , and R1 = 1. Note that Cl is indexed as 0 ≤ l ≤ N and Ri is indexed as 1 ≤ i ≤ ns . The constraint (2.32) is obtained, then, by the minimizer of the functional f = Akl wˆl − u¯k · Akj wˆj − u¯k − λ · (Cl wˆl − Ri u¯i )   (2.34) over the degrees of freedom wˆl with Lagrange multiplier λ. Differentiating with respect to these, we get the equations ∂f ∂f = 2Akl Akj wˆj − 2Akl u¯k − λCl = 0, = Cl wˆl − Ri u¯i = 0 (2.35) ∂ wˆl ∂λ The final linear system to solve, then, that will incorporate these constraints and conserve the mean, is       2Akl Akj −Cl wˆl 2Akl u¯k  · =  (2.36) Cl δlj 0 λ Ri u¯i where δlj is the Kronecker delta which returns one if the two integers l, j are equal, and zero otherwise. This system completely reconstructs the polynomial degrees of freedom while conserving the mean. For clarity, it is noted that the top left block 2Akl Akj is size Np -by-Np , the bottom left Cl δlj is size 1-by-Np , the top right −Cl is size Np -by-1, and the bottom right 0 is a singleton. 28 For the sake of implementation, the left hand side matrix in equation 2.36 will be inverted and stored for each stencil on each element. This is done as a preprocessing stage, once and for all. Then, each reconstruction may be computed by forming the cell average-dependent right hand side and multiplying by the stored matrix. 2.3.4 Weights and the ENO Property Having computed each of the ns reconstructions, the next step is to weigh them in such a way as to preserve high-order accuracy while eliminating non-physical oscillations. The final reconstruction, RW EN O (ξ), then will be ns ns ! X X RW EN O (ξ) = ωi Ri (ξ) = ωi wˆli Ψl (ξ) (2.37) i=1 i=1 where ωi are nonlinear weights which will be determined by considering the smooth- ness of the polynomial reconstructions. Pns To maintain the preservation of the cell average, it is necessary that i=1 = 1, so we the have normalization relationship ω ˜i ωi = Pns (2.38) j=1 ω ˜j The (non-normalized) weights ω ˜ i are of the form λi ω ˜i = . (2.39) ( + σi )r  is a small number close to zero (e.g. 10−8 ), chosen to avoid division by zero. The exponent r is an added nonlinearity to the weights, given to magnify the contribution of smooth reconstructions to the stencil and further penalize reconstructions with oscillations detected. The value r = 4 has been chosen in this work, per suggestion from the original authors. A further, less obvious aspect of the scheme is the stencil biasing parameter λi . 29 This is given as   λ , i=1 1 λi = . (2.40)  1, else The value of λ1 may be adjusted to give an artificial preference to the central stencil, which aids in maintaining high-order accuracy of the reconstruction at the possible risk of allowing oscillations. For the purpose of robustness it would be preferable to have a definite, set value of λ1 , but in practice one finds that for best results it is problem-dependent. An acceptable value is on the order of 103 , while for problems with high-order features (for example, the shock density wave equation in Sec. 3.2.4) that need to be accurately resolved, it may be preferable to increase this to be more in the range of 106 . Conversely, for problems with mostly piecewise linear data such as Sod’s shock tube problem (Sec. 3.2.4) it may be preferable to take a lower value for λ1 , perhaps even as low as 1 (no biasing of the central stencil). In the implemention of the standard WENO scheme [51], the explicit, recon- structed quantities are function values at the cell boundaries. As a result, the smoothness indicators simplify to a non-linear combination of cell averages. Compli- cations arise, then, when considering problems with irregular or unstructured mesh. In the DK reconstruction, however, the entire set of polynomial degrees of freedom w˜i are reconstructed, and so the smoothness indicator may be computed as a quadratic functional of these degrees of freedom, seen in equation (2.41). N Z 2 ∂ |α| X  σi = wi (ξ) dξ (2.41) TE ∂ξ α |α|=1 Here, α = (α1 , . . . , αn ) , |α| = α1 + . . . + αn is a multi-index. Note that because we are working in the reference system, the mesh scaling has been removed from the equation. Eq. (2.41) can be succinctly written as σi = wˆli Σlm wˆm i , l, m = 1, 2, . . . , N (2.42) 30 where the matrix Σlm acts on the degrees of freedom corresponding to the higher order (non-constant) basis functions and is defined as N Z X ∂ |α| ∂ |α| Σlm = Ψl (ξ) · Ψm (ξ)dξ. (2.43) TE ∂ξ α ∂ξ α |α|=1 As all computations are performed over the reference element T E , the matrix Σlm does not depend on the mesh and is therefore a universal smoothness indicator. It may be computed once and for all, using Gaussian quadrature to evaluate the integrals. Thus far, the reconstruction has been presented as for a scalar function u = P j wˆj Pj (x). A discussion in [18], however, details additional features of the scheme when we consider systems of conservation laws. In particular, weighted reconstruc- tions are performed on characteristic variables, in addition to conserved ones in each normal direction to the element, from among which the reconstruction of each com- ponent is chosen. For the results of the one-dimensional Euler system of equations presented in Sec. 3.2.4, this implementation is adopted. Chapter 3 The Hybrid Scheme 3.1 Switching Between Schemes The discontinuous Galerkin scheme offers high-order accuracy on an arbitrary mesh, provided that the underlying solution is sufficiently smooth. Furthermore, the com- putational cost of implementing DG is generally much cheaper than that of a finite volume scheme of similar accuracy, especially if one chooses to use non-linear ENO function reconstructions for the numerical flux. It is desired, then, to maintain the DG scheme as much as possible and to switch to the high-order ENO scheme only when necessary. One such way of accomplishing this is via a limiter. Rather than actually changing the solution, however, we seek to use the limiter as simply a means to judge the smoothness of the solution. The generalized slope limiter presented in [9] provides such a framework. We outline a shock-detecting procedure by first defining the minmod function  1≤i≤n |ai | , |s| = 1  s min Pn 1 m(a1 , ..., an ) = s= n i=1 sign(ai ) (3.1)  0, otherwise, The minmod function will return 0 if the signs of the input ai are not in agreement, 31 32 and otherwise will return the smallest (in magnitude) of the ai . From this, we have the early workings of a possible switch parameter for the Hybrid scheme: given cell averages in a suitably local stencil as inputs to m, we could opt to switch schemes based on the output. The following approach is taken for one dimensional problems: Define left and right interface flux approximations as vlk = u¯kh − m u¯kh − ukl , u¯kh − u¯k−1 k+1 k  h , u ¯ h − u¯ h , (3.2) vrk = u¯kh + m ukr − u¯kh , u¯kh − u¯k−1 ¯k+1  h ,u h − u¯kh . where ukl and ukr refer to the solution at the first (left-most) and last (right-most) gridpoints of the element, respectively. We then make the following judgment over whether a given volume is considered a “DG element” or an “FV element”: If vlk − ukl < 0 and vrk − ukr < 0 , (3.3) then element T k is tentatively considered a DG element. Else, element T k is an FV element. Additionally, once all FV elements have been decided upon, all nearest-neighbors are included as well. While such a choice increases the cost and potentially diminishes the accuracy by decreasing the number of DG elements, numerical observations indicate that it is necessary to include the nearest neighbors to prevent spurious oscillations. From (3.3), then, there arises a single parameter 0 which regulates which scheme is to be used. We will refer to this number as the switching parameter. Generally, it is solution- and equation-dependent, although for the test problems values in the range of 10−2 − 10−4 are sufficient. The method proposed to detect shock regions of the solution is only detailed for one-dimensional problems due to the nature of the limiter involved. For higher dimensions the approach would need to be modified; this is not discussed, as in practice the more difficult task is to find a suitable reconstruction algorithm that 33 can provide high-order, monotonic data on arbitrary mesh. Therefore, any sort of detection approach that errs on the side of over-inclusion of FV regions should initially suffice. 3.1.1 Implementation Algorithm The Hybrid scheme is adapted from the LSERK procedure given in Eq. (1.14) for a fourth-order in time, explicit Runge-Kutta time scheme. We assume that at the outset of the time integration, the solution ukh (x) is known for all elements T k everywhere, regardless of whether T k is in a smooth or non-smooth region. We then simultaneously evolve duDG h = Lh (uh , t) (3.4) dt and uFh V d¯ ¯ h (¯ =L uh , t). (3.5) dt where uDG h is a nodal representation of the solution in the smooth regions, and u¯Fh V are the cell averages of the solution in non-smooth regions. Adapting the same conventions as in (1.14) where additional storage variables are labeled v and inter-stage solutions are labeled p, the scheme is as follows: 1. Before beginning the inner Runge-Kutta loops, compute which elements are in declared shock regions using (3.3) 2. Begin inner Runge-Kutta loops: having all necessary solution data available, compute Runge-Kutta residuals v DG and v¯F V , and the updated solutions pDG and p¯F V 3. Extract cell averages from smooth region (DG) elements 4. Having all necessary cell averages available, compute reconstructions in all shock region (FV) elements 34 5. Repeat (2)-(4) for each of the 5 LSERK stages 6. Return to (1) for the next time step The algorithm will minimize the number of reconstructions needed, as this tends to be the greatest cost. For CFL time stepping restrictions, one may also need to ¯ h operators. consider the eigenvalues of both the Lh and L 3.2 Results In this section we will examine numerical simulations involving the Hybrid scheme. For a variety of transient simulations we will answer the following questions: • Are we able to prevent non-physical oscillations? • Are we able to accurately capture the discontinuity? • What parameter sensitivities do we observe? • How are the results affected by mesh refinement? 3.2.1 Convergence Studies A fundamental property of the reconstruction algorithms is the ability to fully rep- resent smooth functions up to the order of the discretization. To test the O(hN +1 ) convergence properties the reconstructions are applied to the function f (x) = exp [sin (x)] , x ∈ [0, 2π] . (3.6) Convergence studies are performed for each reconstruction and presented in Ta- bles 3.1 and 3.2. These confirm O(hN +1 ) accuracy. 35 Q-ENO DK K L∞ -error Convergence rate L∞ -error Convergence rate 20 1.93e-02 3.07e-02 40 2.22e-03 3.12 4.29e-03 2.84 N = 2 80 2.72e-04 3.03 5.51e-04 2.96 160 3.35e-05 3.02 6.93e-05 2.99 20 6.12e-03 7.58e-03 40 8.15e-04 2.91 4.46e-04 4.09 N = 3 80 6.00e-05 3.76 2.50e-05 4.16 160 3.88e-06 3.95 1.50e-06 4.06 20 7.89e-03 1.51e-02 40 2.48e-04 4.99 6.80e-04 4.47 N = 4 80 5.79e-06 5.41 2.38e-05 4.83 160 1.91e-07 4.92 7.65e-07 4.86 Table 3.1: Convergence study for DK and Q-ENO reconstructions. L∞ error is given for K elements of order N . Q-ENO DK K L1 -error Convergence rate L1 -error Convergence rate 20 1.13e-02 2.58e-02 40 1.42e-03 2.99 3.22e-03 3.00 N = 2 80 1.85e-04 2.95 4.00e-04 3.01 160 2.35e-05 2.97 4.98e-05 3.01 20 6.61e-03 7.60e-03 40 4.32e-04 3.93 3.27e-04 4.53 N = 3 80 2.44e-05 4.15 1.58e-05 4.37 160 1.28e-06 4.25 8.90e-07 4.15 20 4.78e-03 1.67e-02 40 1.19e-04 5.33 6.33e-04 4.72 N = 4 80 4.20e-06 4.82 2.12e-05 4.90 160 1.51e-07 4.80 6.73e-07 4.98 Table 3.2: Convergence study for DK and Q-ENO reconstructions. L1 error is given for K elements of order N . For DK we use a biasing parameter of λ1 = 103 . 36 3.2.2 One-dimensional Advection Equation The one-dimensional advection equation is given as ∂u ∂t + c ∂u ∂x = 0, x ∈ [a, b] (3.7) u(x, 0) = u0 (x) The initial condition u0 is assumed to be periodically extended across [a, b], and so the analytic solution to this equation is u(x, t) = u0 (x − ct). That is, the solution is a traveling wave with speed c. As the solution of the PDE is simply a propagation of the initial condition u0 (x), the smoothness of the solution is completely determined by the initial condition as well. Riemann Problem For our first example, consider the Riemann problem ∂u ∂t + 2π ∂u ∂x = 0, x ∈ [0, 2π]   0, 0 ≤ x < π (3.8) u(x, 0) =  1, π ≤ x ≤ 2π In Figure 3.1 we present results for varying polynomial order and number of elements. The switching criteria 0 = 10−4 is used; there is little sensitivity to this value. It is observed that h-refinement, corresponding to higher number of elements K, appears to be more crucial to resolving the discontinuity. This is not surprising given that the solution in smooth regions is piecewise constant. 37 (N,K) = (1,20) (N,K) = (1,40) (N,K) = (1,80) 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 0 2 4 6 0 2 4 6 (N,K) = (2,20) (N,K) = (2,40) (N,K) = (2,80) 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 0 2 4 6 0 2 4 6 (N,K) = (3,20) (N,K) = (3,40) (N,K) = (3,80) 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 0 2 4 6 0 2 4 6 (N,K) = (4,20) (N,K) = (4,40) (N,K) = (4,80) 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 2 4 6 0 2 4 6 0 2 4 6 Figure 3.1: One-dimensional advection, Riemann problem (see Eq. (3.8)). Terminal time is t = 1. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. Left to right: K = 20, 40, 80. Top to bottom: N = 1, 2, 3, 4. The order of polynomial reconstruction is the same as the order of the DG scheme. A switch tolerance of 0 = 10−4 is used. 38 Advection With Mixed Features Next we consider (3.7) with wave speed c = 1, domain x ∈ [−1, 1] and initial condition  1    6 [G(x, β, z − δ) + 4G(x, β, z) + G(x, β, z + δ)], x ∈ [−0.8, −0.6]   1, x ∈ [−0.4, −0.2]      u(x, 0) = 1 − |10(x − 0.1)|, x ∈ [0, 0.2] (3.9)    1 [F (x, α, a − δ) + 4F (x, α, a) + F (x, α, a + δ)],  x ∈ [0.4, 0.6] 6       0, else where p F (x, α, a) = max(1 − α2 (x − a)2 , 0) (3.10) G(x, β, z) = exp(−β(x − z)2 ) log 2 and z = −0.7, δ = 0.005, β = 36δ 2 ,a = 0.5, and α = 10. The initial condition (3.9), given by Jiang and Shu [30], consists of a Gaussian, a square wave, a triangle and a semi-ellipse. This particular example is very instruc- tive in measuring the Hybrid’s ability both to capture discontinuous data and to accurately resolve smooth features with steep gradients. We perform tests first on a mesh fixed with K = 200 elements, but with varying order of N = 1, 2, 3 and 4. A switch tolerance of 0 = 10−4 is used. Results given in Figure 3.2 indicate that for the low-order N = 1 solution there is excessive smearing of all features. For sufficiently small number of elements K this tends to be the result as the features of the problem are under-resolved. Increasing the polynomial order of approximation, we observe that the solution is well-captured for N = 2 and N = 3. However, spurious oscillations emerge N = 4. This is attributed to an over-ambitious search within the polynomial reconstructions where the underlying data is unsuitable for such a high-order reconstruction. One way to resolve this is to vary the switching parameter 0 = 10−2 . In Figure 3.3 we display results where a looser value, 0 = 10−2 , has been chosen. This in effect allows 39 N=1 N=2 1 1 0.8 0.8 0.6 u(x,t=8) u(x,t=8) 0.6 0.4 0.4 0.2 0.2 0 0 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x N=3 N=4 1 1.2 1 0.8 0.8 0.6 u(x,t=8) u(x,t=8) 0.6 0.4 0.4 0.2 0.2 0 0 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x Figure 3.2: Hybrid solution for one-dimensional advection equation with initial con- dition given by Eq. (3.9). Q-ENO reconstruction is given by O’s; DK is given by X’s. Exact solution is given by solid line. K = 200 elements are used, with varying orders N = 1, 2, 3, 4. Switch tolerance is fixed at 0 = 10−4 . the discontinuous Galerkin method to come back in effect in shock regions, if the shocks have been sufficiently dissipated. While not a perfect solution, it does help to prevent a scenario where inaccurate reconstructions manage to aggravate numerical oscillations, rather than to limit them. 40 N=1 N=2 1 1 0.8 0.8 0.6 0.6 u(x,t=8) u(x,t=8) 0.4 0.4 0.2 0.2 0 0 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x N=3 N=4 1 1 0.8 0.8 0.6 0.6 u(x,t=8) u(x,t=8) 0.4 0.4 0.2 0.2 0 0 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x Figure 3.3: Hybrid solution for one-dimensional advection equation with initial con- dition given by Eq. (3.9). Q-ENO reconstruction is given by O’s; DK is given by X’s. Exact solution is given by solid line. A mesh of K = 200 elements are used, with varying order N = 1, 2, 3, 4. Switch tolerance is 0 = 10−2 41 3.2.3 Burgers’ Equation The one-dimensional Burgers’ equation is given as  2 ∂u u ∂t + 2 = 0, x ∈ [−1, 1] x (3.11) u(x, 0) = u0 (x) Due to the nonlinearity u2 in the flux term of Eq. (3.11), the solution may develop a shock in finite time, even with a smooth initial condition u0 (x). The Hybrid scheme is applied to the one-dimensional Burgers’ Equation with the initial condition u0 (x) = − sin (πx) (3.12) The exact solution to (3.11) with initial condition (3.12) is given by the Lax- Oleinik formula [20]: x − y(x, t) u(x, t) = (3.13) t where for each (x, t) pair, y(x, t) can be most easily computed by solving the for the real zeros y of the function (x − y) − − sin (πy) = 0. (3.14) t Convergence studies are performed by examining the solution at time t = 0.2, which occurs before the shock has formed. For these simulations we use the Hybrid and fix the “shock” regions to be every element in the domain; equivalently, we simply have a finite volume scheme with Q-ENO and DK reconstructions to compute numerical fluxes. These results are presented in Figure 3.3 and 3.4. To better understand the features of the problem, exact and numerical solutions are displayed in Figure 3.4. Error plots for these solutions are given in Figure 3.5. If we were to animate the solution of (3.11)-(3.12), we would observe the positive 42 Q-ENO DK K L∞ -error Convergence rate L∞ -error Convergence rate 20 2.07e-01 1.75e-01 40 4.73e-02 2.13 2.97e-02 2.56 N =1 80 9.15e-03 2.37 9.81e-03 1.60 160 3.51e-03 1.38 4.15e-03 1.24 20 9.70e-02 1.21e-01 40 2.01e-03 2.27 3.49e-02 1.79 N =2 80 3.52e-04 2.51 7.04e-03 2.31 160 4.84e-04 2.86 1.04e-03 2.75 20 1.35e-01 6.82e-02 40 2.10e-02 2.68 1.45e-02 2.23 N =3 80 1.47e-03 3.84 1.48e-03 3.29 160 8.51e-06 4.11 7.85e-05 4.24 20 9.99e-02 8.84e-02 40 1.45e-02 2.78 2.19e-02 4.72 N =4 80 9.95e-04 3.86 2.78e-03 4.90 160 4.70e-05 4.41 1.72e-04 4.98 Table 3.3: Convergence study of Burgers’ equation with smooth initial condition u(x, 0) = − sin (πx). L∞ errors are computed at time t = 0.2. Results are given using a A pure finite volume scheme with Q-ENO and DK reconstructions 43 Q-ENO DK K L1 -error Convergence rate L1 -error Convergence rate 20 6.79e-02 7.08e-02 40 1.67e-02 2.03 1.59e-02 2.16 N =1 80 3.61e-03 2.21 3.51e-03 2.18 160 8.67e-04 2.06 8.14e-04 2.11 20 1.19e-02 1.79e-02 40 1.66e-03 2.85 3.05e-03 2.55 N =2 80 2.81e-04 2.56 4.89e-04 2.64 160 3.80e-05 2.88 6.56e-05 2.90 20 1.38e-02 1.12e-02 40 1.22e-03 3.50 1.57e-03 2.84 N =3 80 8.65e-05 3.82 7.81e-05 4.33 160 4.39e-06 4.30 3.40e-06 4.52 20 9.32e-03 1.91e-02 40 1.01e-03 3.21 2.78e-03 2.78 N =4 80 3.74e-05 4.76 1.85e-04 3.92 160 1.76e-06 4.40 8.19e-06 4.49 Table 3.4: Convergence study of Burgers’ equation with smooth initial condition u(x, 0) = − sin (πx). L1 errors are computed at time t = 0.2. A pure finite volume scheme with Q-ENO and DK reconstructions is used 44 t = 0.2 t=1 1 Q−ENO Q−ENO DK 0.5 DK 0.5 Exact Exact u(x,t) u(x,t) 0 0 −0.5 −0.5 −1 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x Figure 3.4: Solution to Burgers’ equation (3.11)-(3.12) at times t = 0.2 (left) and t = 1 (right). A uniform mesh of order N = 3 and K = 80 elements is used. (left) side of the initial sine wave traveling rightwards, and the negative (right) side traveling leftwards. This forms the shock at x = 0 that is observed in the right graph of Figure 3.4. Naturally, then, the error (in both pre- and post-shock times) is greatest in this central region, as is seen in Figure 3.5. Recall, however, that the principle feature of the Hybrid is to operate on un- structured mesh. The Burgers’ equation (3.11)-(3.12) provides an ideal test problem for this feature, as we know a priori that there will be a shock formed at x = 0. Therefore, we may fix a mesh that refines cell widths around this region. To this end, we prescribe the mapping Vk → sign(V ) · Vk1.5 , k = 1, 2, ..., K + 1 (3.15) where Vk correspond to the vertices of the K elements. Therefore, as we approach the shock region, elements become smaller, providing greater resolution of the solution. 45 t = 0.2 t=1 Q−ENO Q−ENO DK −2 DK 10 −6 10 −4 10 u(x,t) u(x,t) −6 10 −8 10 −8 10 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x Figure 3.5: Error plots corresponding to the simulations of Burgers’ equation in Figure 3.4 are given. t = 0.2 t=1 Q−ENO Q−ENO −4 10 DK −2 DK 10 −4 10 −6 10 u(x,t) u(x,t) −6 10 −8 −8 10 10 −1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1 x x Figure 3.6: Solution to Burgers’ equation (3.11)-(3.12) at times t = 0.2 (left) and t = 1 (right). A non-uniform mesh of order N = 4 and K = 40 elements is used, with element stretching corresponding to Eq. (3.15). 46 3.2.4 One-dimensional Euler Equations The Hybrid scheme is tested on the 1-D system of incompressible Euler equations, given by Qt + F x = 0 (3.16) Q = (ρ, ρu, E)T , F = (ρu, ρu2 + P, (E + P )u)T , with the equation of state given by 1 P = (γ − 1)(E − ρu2 ), (3.17) 2 and γ = 1.4 is taken. For the Q-ENO reconstruction, the density ρ is used to compute the nonlinear smoothness indicator function; see (2.14). Moreover, this reconstruction is applied to the conserved variables. However, the DK method performs reconstruction on characteristic variables. Unless otherwise stated, for all tests we use a switching parameter value of 0 = 10−4 . Sod’s Problem Sod’s shock tube problem presents a simple test of the Hybrid for the one-dimensional Euler Equations. The initial condition is the Riemann problem   (1, 0, 1), x ∈ [0, 0.5) (ρ, u, P )(x, 0) = (3.18) (0.125, 0, 0.1), x ∈ [0.5, 1]. The solution is continuously extended at the x = 0 and x = 1 boundaries, and a terminal time of t = 0.2 is used. An exact solution is available via a Riemann solver. Although computations are done with a switch tolerance of 0 = 10−4 ; experiments reveal that 0 = 10−2 also works at the cost of introducing of small oscillations. In Figure 3.7 numerical solutions of the ρ, u, P and Mach number M are plotted. A fixed mesh of K = 200 elements at N = 1 order is used. The discontinuities are 47 1 0.8 0.8 0.6 0.6 u ρ 0.4 0.4 0.2 0.2 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x 1 0.8 0.8 0.6 0.6 M P 0.4 0.4 0.2 0.2 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x Figure 3.7: Numerical solution to Sod’s problem (3.18). K = 200 elements and N = 1 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. all well-resolved without oscillations. Results of the density ρ for varying order and number of elements are given in Figure 3.8. One observes that increasing both N and K further resolves the discontinuities and avoids oscillations. Lax Problem The Lax problem for the one-dimensional Euler equations (3.16) is given by the initial condition 48 K = 100, N = 1 K = 200, N = 1 1 1 0.8 0.8 ρ(x,t=.2) ρ(x,t=.2) 0.6 0.6 0.4 0.4 0.2 0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x K = 100, N = 2 K = 200, N = 2 1 1 0.8 0.8 ρ(x,t=.2) ρ(x,t=.2) 0.6 0.6 0.4 0.4 0.2 0.2 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x Figure 3.8: Numerical solution to Sod’s problem (3.18). The total number of ele- ments K and polynomial order N is varied. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 49  (0.445, 0.698, 3.528), x ∈ [−5, 0) (ρ, u, P )(x, 0) = (3.19)  (0.5, 0, 0.571), x ∈ [0, 5] The solution is continuously extended at the x = −5 and x = 5 boundaries. A terminal time of t = 1.3 is used. The Lax problem presents a greater challenge than the Sod problem (3.18) due to the shock region (2 < x < 4 at terminal time) initially having zero support. Solutions for a mesh of K = 200 elements at order N = 2 are given in Figure 3.9. An “exact” solution is provided by running the Hybrid with 4000 elements at order 2. In Figure 3.10, a fixed K = 200 elements are used with varying N = 1, 2, 3, 4. The purpose of this simulation is to demonstrate the error introduced by using higher order solutions for the Lax problem. As N increases, the oscillations grow larger. With proper fine tuning of the parameters we may slightly decrease these oscillations, but generally speaking, for this problem it is best to increase K while keeping N low. To demonstrate this, in Figure 3.11 we consider a fixed N = 1 with varying K = 100, 200, 400, 800. It is observed also that the DK scheme is more successful in preventing oscillations than Q-ENO. Shock-Density Wave Problem The Hybrid is applied to the difficult shock-density wave problem, with initial condi- tion (3.20). This problem presents a considerable challenge: as the initial sine waves in the density grow in amplitude, we are at risk smearing the features with excessive dissipation as well as overshooting the correct solution due to aliasing errors.  (3.857143, 2.629369, 10.333333) x ∈ [−5, −4) (ρ, u, P ) (x, 0) = (3.20)  (1 + 0.2 sin(5x), 0, 1) x ∈ [−4, 5] 50 1.5 1.2 1 1 0.8 u ρ 0.5 0.6 0.4 0 −5 0 5 −5 0 5 x x 3.5 3 0.8 2.5 0.6 2 M P 0.4 1.5 0.2 1 0.5 0 −5 0 5 −5 0 5 x x Figure 3.9: Numerical solution to the Lax problem (3.19). K = 200 elements and N = 2 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 51 K = 200, N = 1 K = 200, N = 2 1.2 1.2 1 1 ρ(x,t=1.3) ρ(x,t=1.3) 0.8 0.8 0.6 0.6 0.4 0.4 −5 0 5 −5 0 5 x x K = 200, N = 3 K = 200, N = 4 1.4 1.2 1.2 1 1 ρ(x,t=1.3) ρ(x,t=1.3) 0.8 0.8 0.6 0.6 0.4 0.4 −5 0 5 −5 0 5 x x Figure 3.10: Numerical solution to the Lax problem (3.19). A fixed value of K = 200 grid points with varying N = 1, 2, 3, 4 order is given. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 52 K = 100, N = 1 K = 200, N = 1 1.2 1.2 1 1 ρ(x,t=1.3) ρ(x,t=1.3) 0.8 0.8 0.6 0.6 0.4 0.4 −5 0 5 −5 0 5 x x K = 400, N = 1 K = 800, N = 1 1.2 1.2 1 1 ρ(x,t=1.3) ρ(x,t=1.3) 0.8 0.8 0.6 0.6 0.4 0.4 −5 0 5 −5 0 5 x x Figure 3.11: Numerical solution to the Lax problem (3.19). Solutions of order N = 1 and varying K = 100, 200, 400, 800 are given. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 53 The solution is continuously extended at the x = −5, x = 5 boundaries. A terminal time of t = 1.8 is used. Numerical results for a mesh of K = 400 elements are given in Figure 3.12 at order N = 3 are given in Figure 3.12. An “exact” solution is provided by running the Hybrid with 4000 elements at order 2. Numerical experiments reveal the choice of K = 400 roughly as a minimum for the number of elements needed to best resolve the solution. However, by increasing the switch tolerance to 0 = 10−2 , and permitting the discontinuous Galerkin scheme to more often remain in place, resolution is achieved for lower K and higher N . Examining the density ρ in Figure 3.12, we observe that the Q-ENO reconstruc- tion appears to have slight over- and under-shoots in the naturally oscillatory regions corresponding roughly to 1 < x < 2.5. This behavior grows worse for higher order, as is seen in Figure 3.13. 3.3 Two-Dimensional Observations Thus far the Hybrid scheme has been presented only for one-dimensional problems. Unfortunately, when each scheme was extended to two dimensions, we were unable to robustly retain the ENO property. The main obstacle encountered appeared to be spurious oscillations introduced by demanding high-order polynomial reconstructions in regions where the underlying data was not sufficiently smooth. Nevertheless, we may still verify the ability of each reconstruction to achieve high-order representation of smooth functions. We consider the smooth function u(x, y) = exp [−5 (x2 + y 2 )] , (x, y) ∈ [−1, 1] × [−1, 1] . (3.21) In Figure 3.14 three grids of increasing refinement (K = 146, 568, and 2310 trian- gular elements) are presented. For each mesh we reconstruct Eq. (3.21) from cell averages. In Table 3.5 convergence rates are reported by considering the decay of the L1 error. These rates confirm the desired O(hN +1 ) accuracy for reconstruction 54 2.5 4 2 3 1.5 ρ u 2 1 0.5 1 0 −5 0 5 −5 0 5 x x 1.5 10 8 1 6 P M 4 0.5 2 0 −5 0 5 −5 0 5 x x Figure 3.12: Numerical solution to the Shock density wave problem (3.20). K = 400 elements and N = 3 order are used. Density ρ, velocity u, pressure P , and Mach number M are plotted (one point per element). O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 55 K = 400, N = 2 K = 400, N = 3 4 4 ρ(x,t=1.8) ρ(x,t=1.8) 3 3 2 2 1 1 −5 0 5 −5 0 5 x x K = 400, N = 4 K = 400, N = 5 5 4 4 ρ(x,t=1.8) ρ(x,t=1.8) 3 3 2 2 1 1 −5 0 5 −5 0 5 x x Figure 3.13: Numerical solution to the Shock density wave problem (3.20). K = 400 elements with varying order N = 2, 3, 4, 5 order are used. One point per element is plotted. O’s: Q-ENO reconstruction. X’s: DK reconstruction. Solid line: exact solution. 56 Figure 3.14: Sequence of refined grids with K = 146, 568, 2310 (unstructured) trian- gular elements. Q-ENO DK K L1 -error Convergence rate L1 -error Convergence rate 146 1.39e-01 1.58e-01 N =1 568 3.06e-02 2.23 4.00e-02 2.01 2310 8.64e-03 1.80 6.64e-03 2.56 146 2.03e-02 2.99e-02 N =2 568 2.29e-03 3.21 3.37e-03 3.22 2310 2.61e-04 3.10 2.51e-04 3.70 146 1.30e-02 6.79e-03 N =3 568 6.14e-04 4.50 4.47e-04 4.00 2310 3.74e-05 3.99 2.54e-05 4.09 146 1.04e-02 4.86e-03 N =4 568 2.68e-04 5.39 1.50e-04 5.12 2310 6.72e-06 5.25 4.13e-06 5.12 Table 3.5: Two-dimensional convergence study for DK and Q-ENO reconstructions for the function u(x, y) = exp [−5 (x2 + y 2 )]. L1 error is given for K elements of order N . The schemes attain their desired reconstruction order. √ of smooth functions, where we now take h = 1/ K. 3.4 Concluding Remarks A Hybrid scheme consisting of the discontinuous Galerkin method conjoined to a finite volume method with high-order, ENO-class reconstruction algorithms for the computation of numerical flux, has been presented. A generalized slope limiter is used to detect regions of the solution where shocks may be present. Two reconstruc- 57 tion algorithms, QENO and DK, have been tested in the Hybrid. Unlike traditional ENO and WENO reconstructions, QENO and DK contain no requirement of a struc- tured mesh. Numerical results have been presented for one-dimensional conservation laws which demonstrate the Hybrid’s ability to accurately capture the solution while limiting non-physical oscillations. However, the reconstruction algorithms each failed to maintain monotonicity in the generalization to two-dimensional test problems on unstructured grids. With successful results on nonlinear systems of conservation laws, albeit in the one-dimensional setting, we conclude that the general idea of the proposed Hybrid scheme has much merit, but work remains to be done to extend these successes to higher-dimensional problems. The key feature missing is a suitable reconstruction algorithm. Part II Reactive Infiltration Instability 58 Chapter 4 High-Order Methods for the Reactive Infiltration Instability in an Upwelling Compacting Mantle The generation and segregation of melt from the earth’s interior is a fundamental, but only partially understood, process. Melt is generated by adiabatic decompressional melting during upwelling in the Earth’s upper mantle. Segregation of melt involves two-phase flow in which low viscosity melt percolates through a much more viscous solid matrix [53, 39, 48, 5]. During its upward migration via flow in porous media, melt generated in the deeper part of the upwelling mantle will interact both thermally and chemically with the overlying mantle. The main objective of the present study is to develop high-order numerical methods that can be used to study the chemical interaction between the melt and solid mantle during melt migration in an upwelling mantle column. A fundamental geologic observation that gives insight into the existing theory of melt transport is that the chemical composition of erupted basalt is not in equilibrium with residual mantle at low pressure [40, 57, 19, 35, 34], particularly at diverging plate boundaries. Additionally, it is observed that the mantle is chemically and 59 60 Figure 4.1: Illustration of the formation and coalescence of melt channels. Adiabatic decompressional melting occurs within the mantle. Channels migrate upwards, with impenetrable solid barriers guiding the melt towards the mid-ocean ridge. Picture taken from [31]. lithologically heterogeneous [32]. To preserve the geochemical signature developed at depth, melt must rise from depths of at least 30 km to the surface without extensive re-equilibration with the surrounding mantle. One such mechanism for this involves the reactive infiltration instability (RII) [8, 42, 2], where highly porous regions form due to a feedback between melt flow and dissolution. It has been suggested [42, 15, 33, 2, 1, 56] that high-permeability dunite channels act as conduits for focused flow, where melt may efficiently segregate from its source region while still maintaining its geochemical signature at depth. Numerical studies involving a simple finite difference method have been presented in [56, 55]; however, these results do explicitly not account for the soluble mineral orthopyroxene (opx), whose presence is integral to the formation of dunite channels. Moreover, they do not include an upwelling term and perform tests only for smaller Damk¨ohler numbers. 61 Figure 4.2: The Mid-Atlantic Ridge, a 10,000 km long string of volcanoes. Colors indicate relative age of the crust under the ocean. Red is youngest, with older ages as one approaches the continents. Picture taken from [31]. In this work we utilize the governing equations for reactive dissolution in a de- formable, permeable medium presented in [2, 56], but with the addition of a solid upwelling term and an equation to track the abundance of opx in the matrix. A two-dimensional domain is considered, utilizing the discontinuous Galerkin finite element method in the vertical direction and a Fourier spectral method in the hor- izontal direction. An explicit fourth-order Runge-Kutta time discretization is used. Lateral variations in melt productivity, due to mantle heterogeneity, are represented by perturbations in the melt fraction at the base of the domain. It is observed that the formation and sustaining of melt channels in the simulation box are directly correlated to the perturbation prescribed. In order to circumvent a restrictive time step posed by the reaction term, a formulation of the governing equations is presented where the concentration of the opx mineral is at equilibrium. Numerical solutions tentatively show agreement in 62 the equilibrium and kinetic formulations. A more sophisticated means of tracking of the regions with depleted opx fraction is necessary for accurate treatment of the mass transfer term. 4.1 Governing Equations The governing equations are chosen from [56] and [39], with the addition of an equation to model the presence of the soluble mineral opx. The volume fraction of opx is important as geologists identify melt channel to be areas where opx has completely disappeared. We consider a three phase porous medium consisting of interconnected melt (f) and a matrix comprising of two minerals: olivine (ol) and orthopyroxene (opx). The mass conservation for phase r ∈ S = {f, ol, opx} is given as ∂ (ρr φr ) + ∇ · (ρr φr Vr ) = Γr (4.1) ∂t where φr is the volume fraction, Vr is the phase velocity, and Γr is the bulk solid- to-fluid mass transfer rate of phase r, that is, the reaction rate. All phase fractions sum to unity: φf + φopx + φol = 1 (4.2) We define cr to be the concentration of soluble component in phase r, yielding the solute conservation equation ∂ X X X (ρr cr φr ) + ∇ · (ρr cr φr Vr ) = ∇ · (ρr φr Dr ∇cr ) , (4.3) ∂t r∈S r∈S r∈S where Dr is the diffusion tensor in phase r. The two solid phases ol and opx move with velocity Vol = Vopx = Vs , (4.4) 63 and the fluid velocity is denoted Vf . Conservation of momentum in the fluid is then given by the Darcy’s law which relates P, the pressure in excess of hydrodynamic pressure, to the velocities Vf , Vs . It is given as κφ φf (Vf − Vs ) = − ∇P (4.5) µ where µ is the viscosity of the fluid and κφ is the permeability, expressed as the power law d2 n κφ = φ . (4.6) b f Here, b is a constant and d is the mineral grain size. Both will disappear in the non-dimensionalization. n is an integer that is typically chosen as 2 or 3 for melt migration in the Earth’s mantle; for this work we take n = 3. The system is closed by accounting for conservation of momentum in the solid, given by the equation [39, 54]     ∂P ∂ ∂Vsi ∂Vsj ∂ 2 = η + + ξ − η ∇ · Vs − (1 − φf ) ∆ρgδi3 . (4.7) ∂xi ∂xj ∂xj ∂xi ∂xj 3 Here, η and ξ are the solid phase shear and bulk viscosities, respectively, while ∆ρ is the density difference between the solid and the melt. In order to simplify the governing equations we make the following assumptions 1. A stoichiometric relationship β exists between ol produced and opx consumed by dissolution, so we can reduce the reaction rates Γr to the single degree of freedom Γopx : Γol = −βΓopx (4.8) Γf = − (Γopx + Γol ) = − (1 − β) Γopx (4.9) 2. With a uniformly upwelling mantle, the shear in the matrix is negligible, so we 64 ∂Vsi take ∂xj = 0 if i 6= j. 3. Solid mineral concentrations are constant and uniform, i.e., copx = ceq 0 opx = copx and col = ceq 0 ol = col . 4. Diffusion in the solid is negligible, and isotropic in the fluid, so Dol = Dopx = 0 and Df = D. 5. Porosity is very small (φf  1). 6. The phase densities ρr are constant and everywhere equal except in the buoy- ancy term ∆ρ. We therefore absorb these into the reactive term, Γr → Γr (4.10) ρr Using these assumptions, the mass conservation equations (4.1) become ∂φf + ∇ · (φf Vf ) = − (1 − β) Γopx (4.11) ∂t ∂φopx + ∇ · (φopx Vs ) = Γopx (4.12) ∂t ∂φol + ∇ · (φol Vs ) = −βΓopx (4.13) ∂t whereas the solute conservation equation (4.3) reduces to  ∂ φf cf + φopx c0opx + φol c0ol +∇· φf Vf cf + φopx Vs c0opx + φol Vs c0ol = ∇·(φf D∇cf ) .  ∂t (4.14) We assume a porosity-dependent bulk viscosity 4 ξ0 ξ + η = m, (4.15) 3 φf where ξ0 is a reference viscosity and the parameter m ∈ [0, 1] allows for varying role 65 of the porosity in the bulk viscosity. This assumption, coupled with the simplification of neglecting matrix shear, reduces the momentum equation (4.7) to ! ξ0 ∇P = ∇ ∇ · Vs − (1 − φf ) ∆ρgnz . (4.16) φm f Substituting (4.16) into Darcy’s law (4.5) yields " ! # κφ ξ0 φf (Vf − Vs ) = − ∇ ∇ · Vs − (1 − φf ) ∆ρgnz . (4.17) µ φm f Recalling (4.2), we sum (4.11)-(4.13) and have ∇ · [φf (Vf − Vs )] = −∇ · Vs (4.18) Taking the divergence of (4.17) and equating it to (4.18) then yields the Helmholtz equation φm f p   κφ ∇· (∇p − (1 − φf ) ∆ρgnz ) = , (4.19) µ ξ0 where we have introduced the compaction rate, ξ0 p= ∇ · Vs . (4.20) φm f Darcy’s law (4.17) may be written as κφ φf (Vf − Vs ) = − [∇p − (1 − φf ) ∆ρgnz ] . (4.21) µ Substituting (4.12)-(4.13) in to (4.14), we are left with a simplified equation for the evolution of the soluble material in the melt, ∂ (φf cf ) + ∇ · (cf φf Vf − φf D∇cf ) = − c0opx − βc0ol Γopx .  (4.22) ∂t 66 Collectively, the conservation equations (4.11),(4.12), and (4.22), with the Helmholtz equation (4.19), and Darcy’s law (4.21), form the system of equations for φf , φopx , cf , and p. We need not solve for φol because we solve for φf and φopx (Equation (4.2)). Also, we are not solving for the solid velocity Vs as it can be approximated by a leading uniform upwelling rate. 4.1.1 Reaction Rate At this point, it is necessary to formally define the reaction term Γopx . From there the equations may be non-dimensionalized, so as to be most suitable for numerical simulation. Assuming an unequilibrated fluid concentration cf 6= ceq  f (x) , we take a first-order kinetic approximation Γopx = −αA (φr ) ceq  f (x) − cf . (4.23) The surface area term A (φr ) is a function of the phases φr that moderates how much of the material is available for reaction, with reaction rate α. Introducing the parameter l ∈ [0, 1], we assume A = φlopx . (4.24) The exponent l effectively allows us to control the role, if in any, played by the surface area φopx term. As will be seen later in a linear stability analysis, this term is significant when we examine the convergence of the kinetic to the equilibrium formulations. An alternative formulation given in Sec. 6.3 includes φf in A as well, although the presence of this term does not appear to have a significant effect. The equilibrium concentration is assumed to function of z only over a given length scale L: ∆c ceq 0 f (x) = f (z) = cf + z. (4.25) L 67 Taking a (lower) boundary at z = 0, c0f = cf (z = 0) is the solubility at the base and ∆c is the (constant) change in solubility over length L. The most interesting setting is that of very high mass-fluid transfer rate, i.e., α → ∞. In this limit, the fluid concentration remains at equilibrium (cf = f (z)), with instantaneous melt-solid reaction. An alternative reaction rate for this case, eq Γeq opx , is obtained by substituting cf = cf = f (z) in (4.22) and (4.11). Rearranging terms, we have D∇ · φf ∇ceq − φf Vf · ∇ceq  f f Γeq = (4.26) opx c0opx − βc0ol − (1 − β) ceq f With the assumption of a linearly-varying solubility gradient (4.25), we may simplify the ∇ceq f terms in (4.26). The reaction term may formally be given, therefore, as Γopx = − αφopx (f (z) − cf ) , α < ∞ (4.27) ∂φ D ∂zf − φf w Γeq opx = , α=∞ (4.28) c/∆c − (1 − β) z Lˆ where cˆ = c0opx − βc0ol − (1 − β) ceq f (0) and w = nz · Vf , the vertical component of the fluid velocity. 4.1.2 Boundary Conditions A two-dimensional simulation domain (x, z) ∈ [0, xL ] × [0, zH ] is taken, with the assumption of periodicity of the solution in the x direction. Boundary conditions in z are given as 68 p(0) = 0, φf (0) = φ0f φopx (0) = φ0opx , cf (0) = c0f (nz · Vf ) (0) = w0 , (nz · Vs ) (0) = W0 ∂p ∂cf = 0, =0 ∂z zH ∂z zH (4.29) 4.2 Dimensionless Equilibrium Equations on a Re- active Time Scale For a form more preferable for numerical simulation, we non-dimensionalize the governing equations according to the following (non-dimensionalized variables are denoted with a ∗): (x, z) = (Lx∗ , Lz ∗ ) t = T t∗ φf = φ0f φf ∗ φopx = φ0opx φopx ∗ cf = c0f + ∆ccf ∗ p = P p∗ κφ = κ0 κφ ∗ ˜ s∗ Vs = W0 nz + ∆V V Vf = w0 Vf∗ (4.30) It is noted that the solid velocity Vs has been broken into two components: a uniform ˜ s∗ , where ∆V is assumed to be small. upwelling W0 nz and an additional piece ∆V V 69 δ β η0 φ0f P e Da γ 10−2 0.5 0.133 0.02 ∞ 10 → ∞ 0 → 10−2 2 Table 4.1: Typical values for the dimensionless parameters given in (4.32). The length, pressure and time scales are chosen as s κ0 ξ0 ξ0 ∆V L= m P = m µ φ0f φ0f L n L φ0f d2 T = κ0 = W0 δ b κ0 g∆ρ w0 = ∆V = φ0f w0 δ. (4.31) φ0f µ L is the compaction length [39]. Time adheres to a reactive scale evidenced by the δ in the denominator of T . The system is governed by the dimensionless parameters φ0f φ0f , η0 = φ0opx , δ = (1 − β) ∆c cˆ , cL(φ0opx ) αˆ l (4.32) w0 L W0 Pe = D , Da = φ0f w0 , γ= w0 The γ parameter is most crucial as it effectively controls the rate of upwelling. Typical values for these parameters are given in Table 4.1. Nondimensionalizing the governing equations according to the scales (4.31) yields the dimensionless equations ∂φf γ ∂φf − φ0f Vf · ∇φf + 1 − φ0f φf pφm  =− f − Γopx (4.33) ∂t δ ∂z ∂φopx γ ∂φopx Γopx =− − φ0f Vf · ∇φopx − φ0f φopx φmf p+η (4.34) ∂t δ ∂z 1−β ∂ (φf cf ) 1 δ = −∇ · (φf cf Vf ) + ∇ · (φf ∇cf ) − Γopx (4.35) ∂t Pe ∂  −∇ · δφf ∇p + δpφm  n  1 − φ0f φf φnf   f =− (4.36) ∂z 70 From Darcy’s law the fluid velocity is φf Vf − γφf nz − φ0f Vs = −φnf δ∇p − 1 − φ0f φf nz .     (4.37) To simplify matters, we drop terms of order φ0f and the diffusive term in the fluid concentration equation. The following equations result, from which we will perform numerical simulations. ∂φf γ ∂φf =− + pφm f − Γopx (4.38) ∂t δ ∂z ∂φopx γ ∂φopx Γopx =− +η (4.39) ∂t δ ∂z 1−β ∂ (φf cf ) δ = −∇ · (φf cf Vf ) − Γopx (4.40) ∂t ∂φnf −∇ · δφf ∇p + δpφm  n  f =− (4.41) ∂z φf Vf − γφf nz = −φnf (δ∇p − nz ) (4.42) Denoting the vertical component of fluid velocity w = Vf · nz , we may write the reaction rate as   −Daφlopx (z − cf ) , Da < ∞ Γopx =   . (4.43)  1 1 ∂φf 1−δz P es ∂z − φ f w , Da = ∞ The porosity evolution equation (4.38) was derived by rewriting the Darcy veloc- ity term φf Vf in Eq. (4.11), using (4.21). An alternative evolution equation, (still consistent with those presented) is derived by leaving the Darcy velocity in place, ∂φf δ + ∇ · (φf Vf ) = −δΓopx . (4.44) ∂t We will favor (4.38) over (4.44) for computational purposes, but (4.44) will be helpful later in Sec. 4.3 for computing the one-dimensional steady state. Boundary conditions for (4.38)-(4.43) will be discussed below in Sec. 4.3.1. 71 4.3 One-Dimensional Steady State In absence of an inflow perturbation, the steady-state solution to (4.38)-(4.43) is a function of z only. This steady state will form the basis of a linear stability analysis, as well as a starting point for transient simulations. One frequently occurring term is the volumetric fluid flux (also referred to as the Darcy velocity [4]), q = φf w. (4.45) With this new variable the steady state of the conservation equations (4.44),(4.39), and (4.40) gives rise to the ordinary differential equations dq = − δΓopx . q(0) = q0 (4.46) dz dφopx δ Γopx =η , φopx (0) = 1 (4.47) dz γ1−β dcf δcf − 1 = Γopx , cf (0) = c0f . (4.48) dz q The equations (4.38) and (4.41), which govern the dynamics of compaction, yield the steady state equations dφf γ − δp = − δΓopx , φf (0) = 1 (4.49) dz dφnf  d dp − φnf δ + δp = − , p(0) = p0 . (4.50) dz dz dz The boundary condition on the flux term q is given by (4.42) as   dp q0 = γ − δ − 1 ≈ 1 + γ (4.51) dz Approximate solutions to (4.46)-(4.48) under the assumptions δ, γ  1, are "  # q0 δq02 δηq0 c0f + z q¯ ≈ − 1+ (4.52) 1 − δz Da (1 − δz)4 γ(1 − β)(1 − δz) − δηq0 (c0f + z) 72 " # q0 1 − δc0f δηq0 (c0f − z) c¯f ≈ z − 1− (4.53) Da (1 − δz)2 γ(1 − β)(1 − δz) + δηq0 (c0f − z) δηq0 z 1 δηq02 φ¯opx ≈ 1 − + · δ(1 − β)(1 − δz) Da γ(1 − β)(δz − 1)4 (1 − δc0f ) "  # δηq0 z + c0f (1 − δz) 1−  (4.54) −γ(1 − β)(1 − δz)(1 − δc0f ) + δηq0 z + c0f (1 − δz) The zeroth-order terms in these equations are the solutions to the ODE for the equilibrium formulation, demonstrating convergence as Da → ∞ of the kinetic to equilibrium formulations. Approximate solutions to the compaction equations (4.49)- (4.50) are  1/n δz φ¯f ≈ 1 + q0 (4.55) 1 − δz −nq0 φnf q0 p¯ ≈ 2  ≈− (4.56) (1 − δz) γφf + nφnf (1 − δz)2 4.3.1 Dimensionless Boundary Conditions We would expect the non-dimensional form of the boundary conditions (4.29) to be φf (0) =1 cf (0) = 0 φopx (0) =1 p(0) = 0 (4.57) with Neumann outflow for p and cf . The boundary condition cf (0) = 0 presumes that the concentration is at equilibrium at in flow, while p(0) = 0 effectively prevents flux through the base. These combine to produce a sharp boundary layer near the inflow boundary, evidenced by steep gradients in φf and p. This steepens for increasing Da, although it reaches a critical value in the limit Da → ∞. The presence of such a boundary layer is somewhat artificial. In any event, it is no more arbitrary to assume a scenario absent of the boundary layer than it is with 73 one. Furthermore, the presence of steep gradients in the solution forces us to take higher resolution in the z direction than would otherwise be necessary. Accordingly, the boundary conditions for p and cf may be reformulated so as to essentially remove this boundary layer. These modified boundary conditions computed by evaluating the steady states (4.53) and (4.56) at z = 0: 1+γ cf (0) = − (4.58) Da − δ (1 + γ) 1+γ p(0) = − (4.59) 1 − δz In Figure 4.3 an illustration is given which demonstrates the difference in the base states p¯ and c¯f depending upon which boundary condition is used. 4.3.2 Domain Length The main goal of this project is to study the effect of that dunite channels have on the formation and migration of channels of melt. Dunite channels are marked by the evidence of previously melted opx. Therefore, a relevant domain length to consider is one that extends up to a region where the opx is exhausted. From the steady state (4.54) this is closely approximated ! q0 zH = δ −1 1− γ(1−β) (4.60) q0 + η where as before, q0 = 1 + γ. Typically we take a domain height slightly shorter (∼ 75%) than prescribed in (4.60); when heterogeneity is introduced into the simulation, the height of opx exhaustion will be less than zH due to the formation of melt channels. 74 Figure 4.3: Base states p¯, c¯f for the difference choice of boundary conditions. Top: p(0) = cf (0) = 0. Bottom: modified boundary conditions as given in Eq. (4.58). The domain for cf is shortened to z ∈ [0, 0.02], as outside of this it is essentially parallel to cf = ceq f = z. 75 4.4 Numerical Implementation We perform numerical simulation on the dimensionless equations (4.38)-(4.43). In the vertical direction, we implement a discontinuous Galerkin (DG) finite element method, while in the horizontal direction, we implement a Fourier spectral method. Each of these methods allow for high-order accuracy due to the large basis from which the solutions φf , φopx , cf , p are chosen. In the discussion to follow, we outline the numerical scheme by presenting details of the solution representation, giving examples of the differentiation matrices in each coordinate direction. A fair amount of discussion is given for how the grid points in the DG direction are selected, necessary for the nodal implementation of the different DG operators. We will present the semidiscrete numerical scheme, followed by additional remarks including time integration and filtering for nonlinear stability. 4.4.1 Fourier Spectral Discretization The scheme is discretized in the x-direction by a Fourier spectral method, M X −1 uh (x, ·, t) = u(xi , ·, t)liF S (x), x ∈ [0, xL ] . (4.61) i=0 We assume an even number M of equally spaced horizontal grid points i (4.62) xi = x , M L i = 0, . . . , M − 1, corresponding to the Lagrange interpolating polynomial     1 πN π liF S (x) = sin (x − xi ) cot (x − xj ) . (4.63) M xL xL 76 The (matrix) differentiation operator is given as  h i d FS  (−1)i+j π cot π (xi − xj ) , i 6= j xL xL x (D )ij = l (x) = (4.64) dx j xi  0, i = j. Typically one would compute derivatives with a fast Fourier transform (FFT), but in our case this is insufficient due to the nonlinearity in the Helmholtz equation (4.19). 4.4.2 Discontinuous Galerkin Discretization Much of the background detail of the discontinuous Galerkin method has already been presented in Sec. 1.3. Here we will detail more specifically implementation issues. Building operators for the discontinuous Galerkin method is a bit more complex than the Fourier spectral method, as corresponding matrices are not written in such an easy, closed-form as in (4.64). To build the DG operators we first present a small amount of background material concerning orthogonal polynomials. From this we may compute the grid points, providing a basis for which we can build the nodal- based operators. Jacobi Polynomials (α,β) To begin with, we recall the definition of the classical Jacobi polynomial Pn (ξ) of order n as the solution to the singular Sturm-Liouville eigenvalue problem d d 1 − ξ 2 w(ξ) Pn(α,β) (ξ) + n (n + α + β + 1) w(ξ)Pn(α,β) (ξ) = 0  (4.65) dξ dξ where we are in the (one-dimensional) reference coordinate system ξ ∈ [−1, 1]. The weight function w is given as w(ξ) = (1 − ξ)α (1 − ξ)β . (4.66) 77 (α,β) Multiplying (4.65) through by the polynomial Pm and integrating, we easily ar- rive at the orthogonality of the family of polynomials. We furthermore assume a normalization so that they are orthonormal, Z 1 Pm(α,β) (ξ)Pn(α,β) (ξ)w(ξ)dξ = δmn . (4.67) −1 (0,0) A special case is α = β = 0, resulting in the Legendre polynomials Pn whose weight function is the identity. There is no known way of simply evaluating the Jacobi polynomials, but it may be done via the recurrence relationship (α,β) (α,β) ξPn(α,β) (ξ) = an Pn−1 (ξ) + bn Pn(α,β) (ξ) + an+1 Pn+1 (ξ), (4.68) where the coefficients are given as s 2 n(n + α + β)(n + α)(n + β) an = (4.69) 2n + α + β (2n + α + β − 1)(2n + α + β + 1) α2 − β 2 bn = − . (4.70) (2n + α + β)(2n + α + β + 2) The initial values of the recurrence are available, with Γ (·) as the standard Gamma function: s (α,β) Γ(α + β + 2) P0 (ξ) = 2−α−β−1 (4.71) Γ(α + 1)Γ(β + 1) s (α,β) 1 (α,β) α+β+2 P1 (ξ) = P0 (ξ) ((α + β)ξ + (α − β)). (4.72) 2 (α + 1)(β + 1) Lastly, we mention a useful property of the Jacobi polynomials [58] d (α,β) p (α+1,β+1) Pn (ξ) = n(n + α + β + 1)Pn−1 (ξ). (4.73) dξ The reason we concern ourselves with this discussion of Jacobi Polynomials 78 is found in a well known [16] result connecting Jacobi polynomials and Gaussian quadratures for the approximation of integrals. It can be shown that Z 1 N X +1 f (ξ)w(ξ)dξ = f (ξi )wi (4.74) −1 i=1 (α,β) is exact for polynomials f of degree 2N + 1 if we choose ξi to be the roots of PN +1 (ξ) with corresponding weights wi . To compute ξ and wi we use the recurrence rela- (α,β) tionship (4.68). Setting PN +1 (ξi ) = 0 results in a symmetric tridiagonal eigenvalue problem whose eigenvalues are the Gauss points ξi ; the weights wi may be recovered from the eigenvectors (see [21] for more details). Selecting the DG Nodes All of the machinery presented is utilized in computing the Gauss-Lobatto points, for which we use to represent the solution. For simplicity we present the solution as represented on the reference element ξ ∈ [−1, 1]. A linear transformation then connects ξ to the physical domain z ∈ [0, zH ]. This approach may be expanded to multiple DG elements by simply keeping track of the Jacobians of each (reference to physical) transformation. Although above we have presented a method to evaluate the Gauss points, these points are unsatisfactory as they do not include the end grid points ξ = −1, 1. We prefer to have a nodal representation at these locations in order to more easily im- plement boundary conditions. Instead, for an element of solution order N , we select the N + 1 Gauss-Lobatto nodes. These correspond to the roots of the polynomial  d (0,0) 1 − ξ2 P (ξ) . dξ N Recalling Equation (4.73), we see that the Gauss-Lobatto points are the N − 2 order (α+1,β+1) Gauss points of PN −2 with the end points ξ = −1, 1 added. 79 Gauss−Lobatto Nodes, N = 20 1 0.5 0 ξ −0.5 −1 Figure 4.4: Gauss-Lobatto quadrature nodes for order N = 20. The nodes are given over the one-dimensional reference element ξ ∈ [−1, 1]. The Gauss-Lobatto nodes cluster as O(1/N 2 ) near the boundaries. An illustra- tion of this clustering is given in Figure 4.4 for the value N = 20. The DG Differentiation Matrix Deriving the discontinuous Galerkin derivative operator requires a bit more effort than than was needed for the Fourier differentiation matrix (4.64). In the discussion that follows we compute the local DG differentiation matrix. Unless otherwise stated, all polynomials Pn are assumed to be Legendre polynomials, that is, Pn (ξ) := Pn(0,0) (ξ). Consider the nodal representation of an order N polynomial N X +1 u(ξ) = u(ξ)ljDG (ξ), ξ ∈ [−1, 1] (4.75) j=1 80 with Lagrange interpolating polynomials ljDG based on the N + 1 Gauss-Lobatto points. We seek a matrix Dξ , N +1 ∂u X ∂ljDG = u(ξ) = Dξ u, (4.76) ∂ξ j=1 ∂ξ where the vector u = (u(ξ1 ), ..., u(ξN +1 ))T . The entries of the matrix can be written, then, as ξ  dljDG D ij = . (4.77) dξ ξi To compute the entries of Dξ we recall the Vandermonde matrix introduced in Eq. (1.7), with entries (V)ij = Pj−1 (ξi ). As a consequence of the uniqueness of polynomial interpolation, we have V T l (ξ) = P(ξ), (4.78) where l(ξ) = (l1 (ξ), ..., lN +1 (ξ))T and P = (P0 (ξ), ..., PN (ξ))T . From Eq. (4.78) we differentiate with respect to ξ to get d d VT l(ξ) = P(ξ). (4.79) dξ dξ Combining (4.77) and (4.79), we have T T V T Dξ = Vξ , (4.80) where the matrix V ξ is defined as ξ  dPj V ij = . (4.81) dξ ξi 81 To compute the entries in (4.81) we use the identity (0,0) dPn p (1,1) = n(n + 1)Pn−1 (ξ) (4.82) dξ along with the recurrence relationships from Eq. (4.68). With V ξ and V, we may now compute Dξ from Eq. (4.80), Dξ = V ξ V −1 (4.83) The Dξ matrix presented operates on the reference coordinate (ξ ∈ [−1, 1]), and so the corresponding matrix in physical space (z ∈ [0, zH ]) includes an additional factor from the Jacobian. For a DG element of length h this matrix is ∂z ξ 2 Dz = D = Dξ . (4.84) ∂ξ h The matrix Dz presented above is only part of the story when it comes to the DG differentiation operator. A complete discussion involves the numerical flux. Mass and Stiffness Matrices Later in Sec. 4.4.4 we will present the full implementation of the scheme. There we will use Mξ and S ξ , the (local) mass and stiffness matrices, respectively. These are defined as Z 1 ξ liDG (ξ)ljDG (ξ)dξ  M ij = (4.85) −1 where l are the standard Lagrange interpolating polynomials for the selected set of nodes on which the DG solution is represented, and 1 dljDG (ξ) Z ξ liDG (ξ)  S ij = dξ. (4.86) −1 dξ 82 It can be shown that these matrices are computed easily by the identities −1 Mξ = VV T (4.87) and S ξ = Mξ D ξ . (4.88) As will be seen in the implementation, the stiffness matrix S ξ typically is replaced with the differentiation matrix Dξ , or its physical space representation, Dz . A similar relationship between reference and physical space exists with the mass matrix. For a DG element of length h, we define the mass matrix in physical space, Mz , as h ξ h −1 Mz = M = VV T (4.89) 2 2 We go to the effort of computing the matrix Dξ in reference space, and then transforming to physical space at the end, because in a discontinuous Galerkin dis- cretization the solution is defined piecewise on elements. Therefore we build all DG operators in the setting of the reference element, invoking the transformations as necessary. 4.4.3 Representing the Solution Combining the horizontal and vertical representations of the solution (Eqs. (4.61) and (4.75), respectively), solutions on the k th DG element are represented as M X −1 X N ukh (x, z, t) = ukh (x, z, t)liF S (x)ljk,DG (z) (4.90) i=0 j=0 where liF S , ljk,DG are the Lagrange interpolating polynomials associated with Fourier nodes xi and Gauss-Lobatto nodes zjk , respectively. The computational domain is the (rectangular) tensor product of the Fourier 83  k nodes {xi } and Gauss-Lobatto nodes zj . For implementation purposes, it has been useful to treat the solutions as two-dimensional arrays. In the representation (4.90) this is visualized as             u0,0       u0,1       u0,M −1                         u     u     u   k 1,0 1,1 1,M −1 uh (x, z, t) = u = .. , .. ,··· .. (4.91)       .       .       .                        u     u     u   N,0 N,1 N,M −1 where uj,i = ukh (xi , zjk , t). In such an implementation, operators in z act on columns from the left while operators in x act on rows from the right. For example, derivatives in the x and z directions are discretized as ∂u ∂u ' u (Dx )T , ' Dz u (4.92) ∂x ∂z A global solution (that is, one containing all DG elements) is arranged by vertically concatenating the solution arrays. In this way, global DG operators (e.g. the mass matrix, used in the Helmholtz equation discretization (4.104)) are constructed by considering a sparse, block-diagonal arrangement of each element’s local operator. 4.4.4 Spatial Discretization Prior discussion in Chapter 1 generally outlined the discontinuous Galerkin method for a hyperbolic conservation law, but avoided implementation-specific details. Here we will step through the implementation for the system (4.38)-(4.43). The previously introduced notations for the physical-space mass and differenti- ation matrices is broadened to account for varying DG element sizes. For a given element k of length hk , the notation for these is M(z,k) , D(z,k) 84 so as to denote that they are z-operators applied to the k th discontinuous Galerkin element. We proceed using the strong form of the discontinuous Galerkin method, in- troduced previously in Eq. (1.11). We will evolve the solution on DG element k,   corresponding to physical space z ∈ zBk , zTk . In this case the discrete solution quantities are represented as (φf , φopx , cf , p) ' (φf , φopx , cf , p) . (4.93) In the outline to follow, we attempt to motivate a matrix-vector implementation, invoking the matrices M(z,k) , D(z,k) , and Dx to act on solution vectors (bold-faced quantities; e.g. those in Eq. (4.93)). When consecutive vectors appear, products are computed entry-wise (i.e. on the nodes). Consider the porosity evolution equation (4.38). The local semidiscrete scheme is dφf γ h γ  γ ∗ izTk M(z,k) = − S ξ φf + lk (z) φf − φf dt δ δ δ zBk + M(z,k) pφf m − M(z,k) Γopx (φf , φopx , cf , p), (4.94) where we introduce the notation lk (z) to extract the boundary terms of the solution to compute the numerical flux. In particular, we have lk (zBk ) = (1, 0, . . . , 0)T (4.95) lk (zTk ) = (0, 0, . . . , 1)T . A more preferable form of (4.94) is achieved by multiplying though by the inverse of 85 the mass matrix, yielding dφf γ −1 h k  γ  γ ∗ izTk = − D(z,k) φf + M(z,k) l (z) φf − φf dt δ δ δ zBk + pφf m − Γopx (φf , φopx , cf , p) (4.96) In Eq. (4.96) we have effectively discretized the spatial component of partial differen- tial equation (4.38), and now simply must evolve the solution using a time integration ∗ scheme. We compute the numerical flux term γδ φf with a Lax-Friedrichs flux. The evolution of φopx (Eq. (4.39)) is similar: dφopx γ (z,k)  h k γ γ ∗ izTk (z,k) −1 =− D φopx + M l (z) φopx − φopx dt δ δ δ zBk η + Γopx (φf , φopx , cf , p) (4.97) 1−β Again, a Lax-Friedrichs flux is used. Evolving the fluid concentration requires a bit more effort. Denoting the (local) fluid velocity components Vf = (u, w) , (4.98) we have d (φf cf ) n −1  k zk o = − δ −1 D(z,k) (φf cf w) − M(z,k) l (z) (φf cf w − (φf cf w)∗ ) zTk dt B − δ −1 (φf cf w) (Dx )T − δ −1 Γopx (φf , φopx , cf , p) (4.99) For the numerical flux (φf cf w)∗ we use a global Lax-Friedrichs flux. Numerical experiments have revealed that the highly nonlinear product φf cf w is at significant risk for aliasing errors. We handle this issue with the combination of filtering and higher polynomial order (N ' 20). 86 The only remaining component we need in order to evolve Eqs. (4.96),(4.97), and (4.99) is the compaction rate p. For this we need to solve the elliptic equation (4.41). To accommodate for non-zero Dirichlet boundary conditions p(0) = p0 , (see Eq. (4.58)), we write p˜ p= + p0 (4.100) δ so that p˜(z = 0) = 0, and rewrite (4.41) as ∂φnf   ∂ n dp0 φnf ∇˜ p˜φm −δφm   −∇ · p + f = f p0 +δ φf − := f. (4.101) ∂x dx ∂z Using the previous selections of numerical flux (for the derivatives of φf ), we are able to discretize the right hand side f ' f . From here we write the scheme h   i −M(z,k) φf n p˜ (Dx )T (Dx )T + D(z,k) (φf n q) − pφ ˜ fm z k + lk (z) (q − q ∗ ) zTk = M(z,k) f ,  B (z,k) −1 k zk (z,k) l (z) (p˜ − p˜∗ ) zTk .   q=D u− M (4.102) B For the numerical fluxes p˜∗ and q ∗ we use a central flux and stabilized central flux, respectively. The former choice is made because p has no preferred direction of propagation; the latter is made in order to guarantee invertibility of the operator. These are given as p˜∗ = {{p}} ˜ , q ∗ = {{q}} − τ [p] ˜ (4.103) where {{·}} represents the average of the quantity at the interface and [·] represents the jump. τ is a positive number, which we take to be 1. The role of this parameter is to penalize the solution to disallow jumps in p. ˜ For a greater discussion of the details of the stabilized central flux, as well as other possible choices to use, see [3]. 87 Assembling (4.102) over all k, the resulting linear system is Hp˜ = Mf , (4.104) where M is the global mass matrix, with p˜ and f also representing global arrays. From this point we may utilize a standard linear systems solver. Due to the nonlin- ear φnf permeability terms, H is not symmetric. In the work presented, a matrix-free implementation of the BiCGStab iterative method is used to solve the linear sys- tem. For a preconditioner, we compute the incomplete Cholesky factorization of the operator H with a constant coefficient (φf = 1), and drop tolerance 10−6 . 4.4.5 Time Integration For time integration, a fourth-order, explicit, low-storage Runge-Kutta (LSERK) method is implemented, with CFL number 0.5. As Da tends to infinity, the reaction term Γopx becomes stiff in the time step restriction, perhaps motivating an alternative formulation. However, it has been observed numerically that temporal accuracy is important in resolving the fast-reacting features present in large Da simulations, so to date the high-order Runge-Kutta scheme has been necessary for these cases. The time step restriction in these large (but finite) Da cases motivates us to consider the equilibrium (Da = ∞) setting, which is much easier to integrate in time due to the removal of the Da factor in the reaction term. 4.4.6 Additional Remarks For help in stabilizing the nonlinear terms (particularly those in the fluid velocity Vf ), weak artificial dissipation is introduced to the scheme via an exponential filter in each direction. The relatively weak strength choice s = 16 is used for each direction, with cutoff frequencies of 0 and 6 in the horizontal and vertical directions, respectively. The solution variables φf and p are filtered within each Runge-Kutta 88 stage. The term −1  zk M(z,k) lk (z) zTk B appears often in the numerical flux computations. It is useful to define an operator which evaluates this expression. Lastly, it is noted that the opx concentration φopx is not allowed to attain non- physical, negative values. In the kinetic formulation (Da < ∞) this would, asymp- totically, not occur anyway, as the surface area term in Γopx eliminates dissolution for near-zero values of φopx . In the equilibrium formulation, however, there is no such surface area term and so we set Γopx = 0 when φopx becomes negative, while fixing the opx concentration to 0. In future work we will explore the implementation of a level set to more precisely treat the φopx = 0 region, where dissolution must be turned off. Chapter 5 Validation of Numerical Model Having presented a suitable model for analysis, we next set our sights on validation of the model. In this chapter we detail a linear stability analysis, demonstrate conver- gence of the kinetic to the equilibrium formulation and agreement in the growth rates from the transient solutions with those predicted from the linear stability analysis. 5.1 Linear Stability Analysis for Kinetic Formu- lation We perform a linear stability analysis on the equations (4.38)-(4.43). This study will both help us to better understand the behavior of the system as well as providing a key benchmark that may be used for validation of the transient solution. For each unknown u, u ∈ {φf , φopx , cf , p}, we write ueikx+σt . u(x, z, t) = u¯ + ˜ (5.1) The base state u¯ is assumed to be a function of z only,   1 is a small number, and k ∈ [0, ∞) is the (horizontal) perturbation frequency. σ is the growth rate. With these assumptions we may substitute into the simulation equations (4.38)-(4.43). 89 90 Keeping only O() terms we arrive at the system       ASφf ASφopx AScf ASp φ˜f 1 0 0 0 φ˜f        ˜  0 1 0 0   φ˜opx       0 AOφopx AOcf 0   φopx      = σ     ACφf ACφopx ACcf ACp    c˜f    ¯   0 0 δ φf 0   c˜f         AP φf 0 0 AP p p˜ 0 0 0 0 p˜ (5.2) where the entries are given as ASφf = − γδ D + mφ¯m−1 f p¯ ASφopx = Dalφ¯opx l−1 (z − c¯f ) AScf = −Daφ¯lopx ASp = φ¯m f AOφopx = − γδ D − Dalφ¯l−1 η opx 1−β (z − c ¯f ) η ¯l AOcf = Da 1−β φopx (5.3) = c¯0f nφ¯n−1 0   ACφf f (δ p ¯ − 1) − γ ACφopx = Dalφ¯opx l−1 (z − c¯f ) (1 − δ¯ cf ) = φ¯f (δ p¯0 − 1) − γ φ¯f D − Daφ¯lopx (1 + δ (z − 2¯  n  ACcf cf )) ACp c0f φ¯nf D = δ¯ h   i AP φf = δmφ¯m−1 p ¯ − n (δ p ¯ 0 − 1) ¯n−1 D + φ¯n−1 0 − δ φ¯n−1 p¯00 φ f f f f = δ φ¯f − Dnφ¯0f φ¯fn−1 + (k 2 − D2 ) φ¯nf .  m  AP p As p˜ does not have a growth rate associated with it (p is determined by φf ), the 4-by-4 block eigenvalue problem (5.2) may be simplified by first writing p˜ = − (AP p )−1 AP φf φ˜f (5.4) and then substituting the first and third rows’ p˜-component accordingly. The result- ing 3-by-3 block is then solved easily with a standard eigenvalue solver. The base states φ¯f , φ¯opx , c¯f , p¯ taken as the leading terms of the one-dimensional  91 steady states (4.53)-(4.56),  1/n δz φ¯f = 1 + q0 (5.5) 1 − δz ηq0 δz φ¯opx = 1 − (5.6) γ (1 − β) 1 − δz q0 c¯f = z − (5.7) Da − δq0 q0 p¯ = − (5.8) (1 − δz)2 5.2 Linear Stability Analysis for Equilibrium For- mulation In the event that Da = ∞, the concentration is fixed at equilibrium (cf = z) and the reaction rate has no dependence on φopx . Therefore φf (and in turn, p) are decoupled from φopx in the stability analysis. Accordingly, the linear system (5.2)- (5.3) simplifies to a 2-by-2 block,       A∞ Sφf A∞ Sp φ˜f 1 0 φ˜f    = σ   (5.9) AP φf AP p p˜ 0 0 p˜ with modified fluid concentration entries γ−nφ¯n−1 (δ p¯0 −1) A∞ γ ¯m−1 p¯ − f Sφf = − δ D + mφf δz−1 δ φ¯n (5.10) A∞ Sp = φ¯m f − 1−δz f D. The entries AP φf , AP p for the pressure equation remain the same as the kinetic case. In the same fashion as before, we may reduce (5.9)-(5.10) to a single block by replacing p˜ as a linear combination of φ˜f . 92 5.3 Convergence of Kinetic to Equilibrium For- mulations One measure of validation we may conduct through the linear stability analysis (5.2)- (5.10) is to demonstrate convergence of the kinetic to equilibrium formulations in the limit Da → ∞. This convergence is crucial, as it may justify the usage of the latter formulation for transient numerical simulations. We disregard dependence of the porosity in the viscosity (m = 0) and remove the opx surface area term from reaction (l = 0). Tests are performed with an upwelling rate γ = 0.00025 and increasing Da, presented in Figure 5.1. The top left panel of this figure shows a close match between the growth rates in kinetic and equilibrium formulations for increasing Da. In the top right we measure the convergence of the eigenfunctions φ˜f and p˜, and the eigenvalue σ for a fixed k = 100. That is, we measure the relative differences ˜Da ˜∞ Da φf − φf p˜ − p˜∞ σ − σ∞ Da h, h , ˜∞ |˜ p ∞| σ∞ φf h h where |·|h is the discrete L2 -norm and quantities with a superscript Da refer to kinetic results while those with superscript ∞ refer to the equilibrium. From these results we see a clear convergence of eigenfunctions and eigenvalues between the kinetic and equilibrium cases. For reference, in the bottom left and right panels of Figure 5.1 we plot the eigenfunctions φ˜f , p˜ for the entire spectrum 0 ≤ k ≤ 240 in increments of 20. The eigenfunctions have been normalized so that φ˜f (zH ) = p˜(zH ) = 1. 93 5.4 Matching Transient Growth Rate to Linear Stability Analysis The next validation we perform is to check that the growth of the transient solutions to Eqs. (4.38)-(4.43) match up to the predicted growth rates that were computed in Figure 5.1. In Figure 5.2 these results are presented for the equilibrium formulation. Transient growth rates are approximated by measuring the rate of change of the Fourier coefficient corresponding to the k th mode. We observe a close matching of the transient simulations to the predicted linear growth rates, as shown in the right panel of this figure. One may notice in the transient simulation growth rate results presented in Figure 5.2 that the time scale for equilibrium (top right) and kinetic formulations (bottom right) are wildly different: the former is of order 1 while the latter is on the order of 10−7 . This is due to the kinetic formulation having a time step restriction ∆t ∼ Da−1 . For the results given in Figure 5.2, we have taken Da = 106 . The equilibrium formulation does not carry this restriction. 94 Figure 5.1: Linear stability analysis, as computed in Eqs. (5.2)-(5.10). γ = 0.00025. Top Left: computed growth rates σ vs. perturbation frequency k, for varying Da. Top Right: convergence in φ˜f , p˜, σ as Da → ∞. A fixed perturbation frequency of k = 100 is chosen. Bottom left/right: φ˜f /˜p eigenfunctions for k ∈ [0, 240], plotted in increments of 20. All eigenfunctions have been normalized so that φ˜f (zH ) = p˜(zH ) = 1. 95 Figure 5.2: Comparison of predicted growth rates σ from linear stability analysis (LSA) from Eqs. (5.2)-(5.10) to those actually measured from the transient simu- lations (4.38)-(4.43). Three perturbation frequencies, k = 5, 150, 225 are examined for the transient simulations. Growth rates in the transient simulations are com- puted by considering the rate of change of the magnitude of the Fourier coefficient corresponding to frequency k. Top: Equilibrium formulation. Bottom: Kinetic for- mulation (Da = 106 ) with surface area coefficient l = 0.5. For both cases we take m = 0. Chapter 6 Observations from Numerical Model Having defined the model in Chapter 4 and validated it in Chapter 5, we next turn our attention to exploring features of the model. In this chapter we will first discuss impacts on the stability of the l and m parameters, which control the influence of the opx mineral in reaction and the porosity in bulk viscosity, respectively. Next, we introduce a different non-dimensionalization based on the upwelling velocity. Nu- merical simulations from this dimensionalization demonstrate a clear link between mantle heterogeneity, represented by inflow φf perturbations, and the formation of melt channels. 6.1 Effects of the φm f Viscosity Term Results from the linear stability analysis above demonstrate convergence in the limit Da → ∞ of the kinetic to the equilibrium formulations. Next, we consider the same problem but with varying dependence of the porosity (φm f , m ∈ [0, 1]) in the viscosity term. In Figure 6.1 the predicted growth rates from the linear stability analysis are given 96 97 Figure 6.1: Growth rates σ vs. horizontal perturbation frequency k from linear stability analysis. γ = 0.00025. Left: fixing l = 0, we vary the porosity component φm 6 f in bulk viscosity. A reaction rate constant of Da = 10 is used, which is essentially converged to the equilibrium case. Right: fixing m = 0, we vary the surface area component φlopx in reaction rate. A reaction rate constant of Da = 104 is used. Taking both l, m → 1 stabilizes the system. with varying m, and fixed γ = 0.00025, Da = 106 , and l = 0, the exponent in the opx mineral in the surface term. It is observed that the growth rate σ varies linearly with the m exponent, so that including the presence of porosity in the viscosity term acts to stabilize the system. 6.2 Effects of the φlopx Surface Area Term In the transient numerical simulations presented in Sec. 6.4 it was discovered that the presence of a φopx -dependent surface area term in the kinetic formulation was crucial in acting to turn off reaction once the mineral opx had been dissolved. For the equilibrium formulation this term was no longer present, so that in order to match up the transient solutions we had to include a condition to manually turn off reaction once opx is exhausted (Eqs. (6.19) and (6.20)). 98 Extrapolating these results to our rescaled equations, it is significant to observe that the convergence of kinetic to equilibrium formulations in the linear stability analysis, shown in Figure 5.1, occurs with opx exponent l = 0, i.e., no inclusion of the φopx surface area term in the reaction rate. In the right panel of Figure 6.1 we plot the growth rate for varying l ∈ [0, 1] and reaction rate constant Da = 104 . As with the φf term in the viscosity, the φlopx surface area term stabilizes the equation, with increasing effectiveness as l → 1. This feature is quite minimal, however, and for sufficiently high Da the dispersion relationships for varying l are nearly indistinguishable. 6.3 Dimensionless Governing Equations on an Up- welling Time Scale Previous numerical results have been given with a reactive time scale, specified in Eq. (4.31). In this section we present a slightly different non-dimensionalization, where the time scale is normalized to the upwelling velocity. The motivation behind this dimensionalization is to explore the effects of a source perturbation, simulat- ing mantle heterogeneity, on the sustained formation of melt channels. Transient numerical results will follow in subsequent sections. The surface area and viscosity exponents (l and m, respectively) have been fixed at 1. We now assume length, pressure and time scales as q κ0 ξ0 L = µφ0f ξ0 P = T (6.1) L T = W0 n 2 (φ0f ) d κ0 = b 99 δ β η0 φ0f P e Da R 10−2 0.5 0.133 0.02 ∞ 2 5 10 → 10 , ∞ 100 Table 6.1: Typical values for the dimensionless parameters given in (6.2). The system is described by the following dimensionless parameters φ0f φ0f , η0 = φ0opx , δ = (1 − β) ∆c cˆ , αˆcLφ0 (6.2) W0 L κ0 g∆ρ P es = D , Das = DaR = W0opx , R= µφ0f W0 where cˆ = c0opx − βc0ol − (1 − β) c0f . The parameter δ effectively controls the slope of the solubility gradient (f (z) in (4.25)), with stoichiometric dissolution constant β. The reaction rate constant is given by the product Das = DaR, where Da is a measure of the dissolution rate relative to the rate of melt flow in the porous rock and is referred to as the Damk¨ohler number. The Peclet number P es is a measure of the relative rate of solid upwelling to the rate of diffusion in the melt. In these simulations, we neglect diffusion and so take P es = ∞. Finally, the parameter R is a measure of the ratio of the fluid velocity to solid upwelling velocity. Therefore, higher R corresponds to less upwelling. Typical values which we choose for these parameters are given in Table 6.1. Dropping the ∗ superscripts on dimensionless variables, from these simplifications we arrive at the following dimensionless set of equations for φf , φopx , cf , p, and Vf : ∂φf ∂φf ∆V ˜ = −φ0f φf 2 p − − Vs · ∇φf + φf p − δΓopx (6.3) ∂t ∂z W0 ∂φopx ∂φopx ∆V ˜ η0 δ =− − Vs · ∇φopx − φ0f φf φopx p + Γopx (6.4) ∂t ∂z W0 1−β   ∂φf cf φf = −R∇ · cf φf Vf − ∇cf − Γopx (6.5) ∂t P es ∂  n −∇ · φnf ∇p + φf p = −R φf 1 − φ0f φf    (6.6) ∂z 100 ˜s ∆V V − φnf ∇p − R 1 − φ0f φf nz   φf Vf = φf nz + (6.7) W0 To further simplify matters we disregard terms of order φ0f and ∆V . Doing so yields ∂φf ∂φf =− + φf p − δΓopx (6.8) ∂t ∂z ∂φopx ∂φopx η0 δ =− + Γopx (6.9) ∂t ∂z 1−β   ∂φf cf φf = −R∇ · cf φf Vf − ∇cf − Γopx (6.10) ∂t P es ∂φnf −∇ · φnf ∇p + φf p = −R   (6.11) ∂z φf Vf = φf nz − φnf (∇p − Rnz ) (6.12) The reaction term is   −Das φf φopx (z − cf ) , Das < ∞ Γopx =   (6.13)  1 1 ∂φf 1−δz − Rφf w , Das = ∞ P es ∂z The porosity φf has now been included in the surface area term in reaction; this component does not significantly affect numerical results. The equations (6.8) - (6.13) form the system that we will explore in depth with numerical simulation. As in the previous non-dimensionalization, we modify the cf and p boundary conditions to remove the artificial compaction layer near inflow. These are given as 1+R cf (0) = (6.14) δ (1 + R) − Das   δ 1+R + Da c (0) , Da < ∞ 1+nR s f s p(0) = (6.15)  δ 1+R − 1 − R ,  Das = ∞ 1+nR As before, φopx (0) = 1 and φf (0) = 1 + h(x, t), for a small, prescribed perturbation h. 101 The maximum domain length of interest zH (corresponding to the approximate location of φopx = 0 in the one-dimensional steady state) is 1−β zH = δ −1 (6.16) (1 + R) η + 1 − β 6.4 Transient Numerical Results for the Kinetic Formulation In this section we present numerical simulations on the dimensionless equations (6.8)-(6.13). A strong correlation between mantle heterogeneity, represented by a perturbation in φf at inlet, and the formation of melt channels is discovered. We investigate this and other behavior for the kinetic formulation. In the subsequent section, results are linked with the equilibrium formulation. 6.4.1 Correlation of Inflow φf - Perturbations to Outflow Melt Channels One remarkable feature observed is a strong correlation of the mantle heterogeneity (represented as perturbations in the lower φf boundary condition) to the formation of channels higher in the domain. To demonstrate this we consider a transient simulation with the following prescribed boundary conditions: φf (x, z = 0, t) =     1 + Noise(5), 0 ≤ t ≤ 7.33  1 + Noise(10), 7.33 < t ≤ 11    1 + 0.05 exp −100 (x − 1/6)2 + S(t) (x − 1/2)2 + (x − 5/6)2 , 11 < t ≤ 20    (6.17) 102 where Noise(k) is Gaussian white noise in the first k Fourier frequencies with variance 10−3 , and the switch function S(t) returns 0 if the integer part of t is odd, and 1 otherwise. In Figure 6.2, results are shown. In the bottom panel, the perturbed incoming φf , corresponding to (6.17), is shown, while in the top panel, the outgoing φf is given. We may observe a direct correlation of the two, although it is off by a few (roughly 2) dimensionless time units. This is due to the time it takes for heterogeneities at the base to rise to the top of the simulation domain. Initially, the smoother perturbation of the system (0 ≤ t ≤ 7.33) leads to well- defined, but time-varying channels at outflow. For the higher frequency noise on the time interval 7.33 < t < 11, narrower, more defined channels are observed at outflow. Finally, the alternating 2-3 Gaussian bumps, prescribed on the time interval 11 < t ≤ 20, are observed directly on the outflow. As the middle bump disappears, the middle melt channel at outflow weakens but does not altogether disappear. To better understand the transient behavior of this system with source perturba- tion as prescribed in Eq. (6.17), we plot the solution variables φf , φopx , p at the early time t = 2.8, which corresponds to the first perturbation in Eq. (6.17). For this study, however, we consider the much faster reacting rate of Da = 1000. These plots are given in Figure 6.3. For the solution time chosen, the φf inflow perturbation contains Gaussian white noise in the first five Fourier modes. It is observed that there are five channels appearing in the top of the simulation box, confirming the previous observations of relation between channel formation to source perturbation. Although the channel inflow has porosity around one, at the outflow the channels of melt have formed of porosity nearly 1.4, with accompanying compaction regions where the porosity is less than 0.5. The middle picture shows the opx fraction re- maining, which corresponds as expected with the porosity (less fluid, more solid, and vice versa). The rightmost picture shows the compaction rate. 103 Figure 6.2: Comparison of input (bottom) perturbation to output (top) channels formed. Porosity φf is shown for a reaction rate of Da = 10 up to a dimensionless time of t = 20. Porosity perturbations from Eq. (6.17) are used: For the first 7.33 time units, a low-frequency perturbation with noise in the first five Fourier modes is used. For 7.33 < t < 11, a higher frequency perturbation, with noise in the first 10 Fourier modes, is given. For 11 < t < 20, a regulated series of three Gaussian bumps is prescribed, with the middle disappearing every other time unit. 104 Figure 6.3: Transient simulation with perturbed φf inflow as pictured in the bottom panel of Figure 6.2. For this case, however, we take Da = 1000. The solution is examined at time t = 2.8, which corresponds to the perturbation in the first time regime in Eq. (6.17). Five channels of melt are observed, corresponding to the perturbation noise at inflow contained with the first five Fourier modes. Melt channels in φf correspond one-to-one with dissolved regions of φopx . 6.4.2 Channels Deepen as Da → ∞ An important verification for the numerical solution is the effects of the reaction constant Das = DaR, where Da is the Damk¨ohler number and the constant R effectively is the ratio of the fluid velocity to the solid (upwelling) velocity, which we fix at R = 102 . Increasing Da causes faster dissolution of the mineral opx, which should lead to deeper melt channels. Recalling the transient simulation with varying perturbation presented in Figures 6.2 and 6.3, we now instead consider the effects of increasing Da. From Figure 6.3 it is observed that the most interesting behavior occurs near the top of the simulation box, where melt channels have formed. Next, we consider the solution at outflow at the same time as in Figure 6.3, but with varying Da. These results are presented in Figure 6.4. We observe that for higher Da, there exist deeper channels of melt and greater dissolution of opx. For the most highly reacting cases of Da = 500, 1000, we 105 observe small Gibbs oscillations in the opx fraction. These may possibly be induced by the steep gradients (in the x direction) of φopx , and therefore, Γopx . Moreover, in the vertical direction, our choice of fixing φopx to 0 in the event that it go slightly below this threshold introduces a slight kink in the solution, leading to oscillations when evaluating the upwelling derivative. The results given in Figures 6.2, 6.3, and 6.4 present a convincing case as to the correlation between channel formation and mantle heterogeneity, as well as the effects of increasing the Damk¨ohler number leading to more pronounced channel formations. This latter observation is important because we mainly are interested in the case of Da → ∞, with the hypothesis that high-permeability melt regions may act as conduits for melt migration to the surface. To better analyze this hypothesis, we next will consider the scenario of a fixed perturbation in φf at inflow, yielding a single melt channel. 6.4.3 Observations From a Single Melt Channel A single channel of melt can be induced with the constant φf boundary perturbation φf (x, z = 0, t) = 1 + 0.1 exp −100 (x − .5)2 .   (6.18) With a constant perturbation, the system will ultimately reach a steady state. For transient numerical simulations, this state is typically reached at a dimensionless time in the neighborhood of the length of the z domain. In Figure 6.5, a transient solution for Da = 500 is given at dimensionless time t = 2. Here, the propagation of the boundary perturbation is observed, with compaction channels forming on the outside of the center channel of melt. A domain (x, z) ∈ [0, 1] × [0, 2] is used. A few interesting features are observed near the top of the computational domain. As the opx is exhausted, the reaction rate tends to 0 as there is no mineral left to dissolve. This in turn contributes to a slight drop-off in the porosity, from which 106 φ at outflow, varying Da f 1.4 1.3 1.2 1.1 1 φopx 0.9 0.8 0.7 Da=10 Da=50 0.6 Da=100 Da=500 0.5 Da=1000 0.4 0 0.2 0.4 0.6 0.8 1 x φ at outflow, varying Da opx 0.35 Da=10 Da=50 0.3 Da=100 Da=500 Da=1000 0.25 0.2 φopx 0.15 0.1 0.05 0 0 0.2 0.4 0.6 0.8 1 x Figure 6.4: Solution φf and φopx at outflow (z = 3.66) for varying Da. The same time t = 2.8 as in Figure 6.3 is chosen, with perturbation defined in as in Eq. (6.17) and plotted in Figure 6.2. Higher Da produces deeper channels of melt and opx dissolution. 107 the compaction rate also rises. In the middle and lower panels of Figure 6.5 we plot focus on vertical profile centered at the channel midpoint x = 0.5 to better visualize this behavior. As the profiles in Figure 6.5 demonstrate a significant change in the solution behavior when φopx tends to 0, it begs the question of what happens further out (i.e. higher in the z domain). A very interesting regime arises in this setting, as is seen in Figure 6.6. Here, we observe that the melt channel centered at x = 0.5 actually bifurcates in a “Y”-shape form soon after opx exhaustion. Moreover, there appear to be secondary channels forming near the top of the z domain, in the vicinity of x = 0.25 and x = 0.75. We believe this behavior can be explained by the relationship between the opx fraction φopx and reaction rate Γopx . As in the shorter z domain case presented in Figure 6.5, the depletion of opx reduces the reaction rate, diminishing the fluid production and therefore φf as well. This occurs around z = 2. As we move upwards then, the regions of highest porosity are only able to sustain themselves by feeding on the opx that remains outside the middle of the now-dissolved channel. Therefore, the original melt channel splits for these regions. Meanwhile, the compaction regions (areas of low fluid fraction) on the outside of the original channel induce secondary channels as fluid from the lower half of the domain flows outside and around these compaction regions. To further understand the melt channel bifurcations, we present in Figure 6.7 the fluid velocity Vf = (vx , vz ) fields and corresponding streamlines. 6.5 Transient Numerical Results for the Equilib- rium Formulation The results presented thus far have been for the kinetic formulation, i.e., finite re- action constant Da. The most physically relevant (and interesting, in lieu of the results given) setting is in the limit Da → ∞. For large but finite Da, however, the 108 φf(x=0.5) φopx(x=0.5) 1.6 1 1.4 0.5 1.2 1 0 0 1 2 0 1 2 z z p(x=0.5) Γopx(x=0.5) 0 −0.8 −1 −100 −1.2 −200 −1.4 −1.6 −300 0 1 2 0 1 2 z z Figure 6.5: Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Top row: Two-dimensional profiles for Da = 500. Middle and bottom rows: solutions at in the channel midpoint x = 0.5 as a function of z. Once the opx mineral is exhausted, the porosity drops and compaction rate rises sharply. 109 Figure 6.6: Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Reaction constant is Da = 500. Soon after the opx mineral is exhausted, the melt channel bifurcates into a “Y” formation. time step in the explicit Runge-Kutta time integration will scale as ∆t ∼ (Da)−1 . To date, implementing alternative time stepping schemes (such as Strang splitting with an implicit method for the reaction terms) have not been fruitful, as temporal accuracy appears crucial in the Γopx term. However, we recall from earlier that it is possible to formulate an equilibrium (Da = ∞) version of the equations, with fixed cf (x, z, t) = z. The reaction term for this case is   1 1 ∂φf Γopx = 1−δz P es ∂z − Rφf w , Da = ∞, (6.19) which no longer carries a stiff time step restriction. Here, w = Vf · nz refers to the vertical component of fluid velocity. However, this formulation is incomplete. Recall that the kinetic reaction rate has a surface area term φopx , which smoothly reduces the reaction to 0 as the opx mineral is exhausted. The equilibrium formulation no longer contains this factor, however, which leads to continual melt production despite 110 Fluid velocity streamlines 4 3.5 3 2.5 2 z 1.5 1 0.5 0 0 0.2 0.4 0.6 0.8 1 x Figure 6.7: Top: fluid velocity vector Vf = (vx , vz ), for profiles presented in Figure 6.6. For vx , negative (positive) values represent leftwards (rightwards) flow. Before the channel bifurcates, the fluid tends to flow more inwards towards the channel, while after the bifurcation it flows outwards. The vertical component vz resembles the porosity. Bottom: fluid velocity streamlines. The fluid tends to avoid low- porosity compaction regions. 111 there not being a dissolvable source (opx) remaining. One fairly straightforward, if not naive, way to deal with this is to simply take Γopx (x, z, t) = 0, if φopx (x, z, t) = 0, (6.20) and accepting Eq. (6.19) otherwise. As before we set any negative values of φopx to 0. In Figure 6.8, results are presented for the same problem as in Figure 6.6, but with Da = ∞ and Γopx modified as described in Eqs. (6.19)-(6.20). We observe that the two solution sets appear very similar to another, although not surprisingly, the equilibrium case contains deeper melt channels. This is not to suggest, however, that such an approach is completely foolproof. By simply setting the reaction term to 0 on gridpoints where the opx is exhausted we have made an approximation to the true location of the channel interface, which currently prevents convergence of the kinetic to equilibrium versions of the code. In Figure 6.9 a comparison of varying Da values is explored, demonstrating that although the equilibrium solution closely resembles the kinetic, they are not in com- plete agreement. In future work we will explore more sophisticated models of tracking the φopx = 0 regions in pursuit reconciling the two formulations. 6.6 Concluding Remarks In this work, a high-order numerical scheme consisting of a discontinuous Galerkin finite element method and Fourier spectral method have been used to model melt generation via the reactive infiltration stability in an upwelling compacting man- tle. The numerical discretization has been verified by matching growth rates pre- dicted from a linear stability analysis to transient behavior of the scheme. Two non-dimensionalizations, one based upon a reactive time scale and one based upon an upwelling time scale, have been considered. A strong correlation between mantle heterogeneity and channel formation has 112 Figure 6.8: Transient solution with fixed boundary perturbation given in Eq. (6.18). Terminal time is t = 2. Reaction constant is Da = ∞. Numerical results are very similar to Figure 6.6, but are available at considerably less computational expense. been discovered, aiding in the of magma transport in the Earth’s mantle. Unlike previous models, we account for the dissolvable mineral opx present, and incorporate it as a surface area term in the reaction rate. From this we observe that channels bifurcate in the event of exhaustion of opx mineral, a previously unseen phenomenon. In the limit of the reaction rate constant relative to melt velocity going to infinity (that is, the equilibrium formulation), we can reformulate the reaction term in a fashion that eliminates the opx-dependent surface area, effectively decoupling the system into the two variables φf and p. Still, the opx fraction plays a crucial role in determining where reaction must be turned off. Future areas of research should include a more precise treatment for the reaction rate in the equilibrium formulation, so as to ensure that reaction is not occurring when opx has been depleted. Finally, the kinetic formulation is hindered by a restrictive time step that scales as the order of the reaction constant. Reconciling this issue, perhaps with a type of 113 φf(x,z=3.6684) φopx(x,z=3.6684) 1.6 0.5 Da=100 Da=100 Da=250 Da=250 1.4 Da=500 Da=500 0.4 Da=∞ Da=∞ 1.2 0.3 1 φopx φf 0.8 0.2 0.6 0.1 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x Γopx(x,z=3.6684) p(x,z=3.6684) 50 0 Da=100 Da=100 Da=250 Da=250 0 −0.2 Da=500 Da=500 Da=∞ Da=∞ −50 −0.4 −100 −0.6 Γopx p −150 −0.8 −200 −1 −250 −1.2 −300 −1.4 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 x x Figure 6.9: Comparison of varying Da behaviors at solution outflow z ≈ 3.66, time t = 2, for transient solution with fixed boundary perturbation given in Eq. (6.18). The equilibrium solution closely resembles the kinetic. 114 semi-implicit time stepping scheme, is preferable for lowering computational expense. This improvement, coupled with a more careful treatment of opx-free regions in the equilibrium formulation mentioned above, will enable us to better demonstrate convergence of the two formulations as the reaction constant approaches infinity. Bibliography [1] E. Aharonov, M. Spiegelman, and P. Kelemen. Three-dimensional flow and reaction in porous media: Implications for the earth’s mantle and sedimentary basins. J. Geophys. Res., 102:821–834, 1997. [2] E. Aharonov, J. Whitehead, P. B. Kelemen, and M. Spiegelman. Channeling instability of upwelling melt in the mantle. J. Geophys. Res., 100:433–450, 1995. [3] D. Arnold, F. Brezzi, B. Cockburn, and L. Marini. Unified analysis of dis- continuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal., 39:1749–1779, 2002. [4] J. Bear. Dynamics of Fluids in Porous Media. American Elsevier Pub. Co., 1972. [5] D. Bercovici, Y. Ricard, and G. Schubert. A two-phase model for compaction and damage 1. General theory. J. Geophys. Res., 106:8887–8906, 2001. [6] J. Butcher. The Numerical Analysis of Ordinary Differential Equations: Runge- Kutta and General Linear Methods. John Wiley & Sons, 2003. [7] M. H. Carpenter and C. A. Kennedy. Fourth-order 2n-storage Runge-Kutta schemes. Nasa tm 109112, NASA Langley Research Center, 1994. [8] J. Chadam, D. Hoff, E. Merino, P. Ortoleva, and A. Sen. Reactive infiltration instability. J. Appl. Math, 36:207–221, 1986. 115 116 [9] B. Cockburn. Discontinuous Galerkin methods for convection-dominated prob- lems. In High-Order Methods for Computational Physics. 1999. [10] B. Cockburn, S. Hou, and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws iv: The multidimensional case. Math. Comp., 54:545–581, 1990. [11] B. Cockburn, S. Lin, and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws iii: One- dimensional systems. J. Comput. Phys., 84:90–113, 1989. [12] B. Cockburn and C.-W. Shu. The Runge-Kutta local projection p1 - discontinuous Galerkin method for scalar conservation laws. RAIRO Mod´el. Math. Anal. Num´er., 25:337–361, 1991. [13] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws ii: General framework. Math. Comp., 52:411–435, 1998. [14] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinu- ous Galerkin finite element method for conservation laws v: Multidimensional systems. J. Comput. Phys., 141:199–224, 1998. [15] M. J. Daines and D. L. Kohlstedt. The transition from porous to channelized flow due to melt/rock reaction during melt migration. Geophys. Res. Lett., 21:145–148, 1994. [16] P. Davis and P. Rabinowitz. Methods of Numerical Integration. Academic Press, 1975. [17] M. Dumbser and M. Kaser. Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic problems. J. Comput. Phys., 221:693–723, 2007. 117 [18] M. Dumbser, M. Kaser, V. A. Titarev, and E. F. Toro. Quadrature-free non- oscillatory finite volume schemes on unstructured meshes for nonlinear hyper- bolic systems. J. Comput. Phys., 226:204–243, 2007. [19] D. Elthon and C. M. Scarfe. High-pressure phrase equilibria of a high-magnesia basalt and the genesis of primary oceanic basalts. Am. Mineral., 69:1–15, 1984. [20] L. C. Evans. Partial Differential Equations, pages 205–206. American Mathe- matical Society, 1998. [21] G. Golub and J. Welsch. Calculation of gauss quadrature rules. Math. Comp., 23:221–230, 1969. [22] D. Gottlieb and J. Hesthaven. Spectral methods for hyperbolic problems. J. Comp. Appl. Math., 128:83–131, 2001. [23] E. Hairer, S. N. rsett, and G. Wanner. Solving Ordinary Differential Equations I. Nonstiff Problems. Springer Series in Comput. Mathematics. Springer-Verlag, 1987. [24] E. Hairer and G. Wanner. Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer Series in Comput. Mathematics. Springer-Verlag, 1991. [25] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy. Uniformly high order essentially non-oscillatory schemes, iii. J. Comput. Phys., 71:231–303, 1987. [26] J. Hesthaven. From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex. SIAM J. Numer. Anal., 35:655–676, 1998. [27] J. Hesthaven, S. Gottlieb, and D. Gottlieb. Spectral Methods for Time- Dependent Problems. Cambridge University Press, 2006. 118 [28] J. Hesthaven and T. Warburton. Nodal Discontinuous Galerkin Methods: Al- gorithms, Analysis, and Applications. Springer, 2008. [29] G. Jiang and C.-W. Shu. On a cell entropy inequality for discontinuous Galerkin methods. Math. Comp., 62:531–538, 1994. [30] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted eno schemes. J. Comput. Phys., 126:202–228, 1996. [31] P. B. Kelemen. The origin of the land under the sea. Sci. Am., pages 52–57, Feb. 2009. [32] P. B. Kelemen, G. Hirth, N. Shimizu, M. Spiegelman, and H. Dick. A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges. Phil. Trans. R. Soc. Lond., 355:282–318, 1997. [33] P. B. Kelemen, J. A. Whitehead, E. Aharonov, and K. A. Jordahl. Experiments on flow focusing in soluble porous media, with applications to melt extraction from the mantle. J. Geophys. Res., 100:475–496, 1995. [34] R. J. Kinzler and T. L. Grove. Primary magmas of mid-ocean ridge basalts. 2 applications. J. Geophys. Res., 97:6907–6926, 1992. [35] E. Klein and C. H. Langmuir. Global corrections of ocean ridge basalt chemistry with axial depth and crustal thickness. J. Geophys. Res., 92:8089–8115, 1987. [36] R. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge Uni- versity Press, 2002. [37] X.-D. Liu, S. osher, and T. Chan. Weighted essentially non-oscillatory schemes. J. Comput. Phys., 115:200–212, 1994. [38] H. Luo, J. D. Baum, and R. Lohner. A Hermite WENO-based limiter for discon- tinuous Galerkin method on unstructured grids. J. Comput. Phys., 225:686–713, 2007. 119 [39] D. McKenzie. The generation and compaction of partially molten rock. J. Petrol., 25:713–765, 1984. [40] M. J. O’Hara. Primary magmas and the origin of basalts. Scot. J. Geol., 1:19–40, 1965. [41] C. F. Ollivier-Gooch. Quasi-eno schemes for unstructured meshes based on un- limited data-dependent least-squares reconstruction. J. Comput. Phys., 133:6– 17, 1997. [42] P. Ortoleva, E. Merino, C. Moore, and J. Chadam. Geochemical self- organization, i, reaction-transport feedbacks. Am. J. sci., 287:979–1007, 1987. [43] J. Qiu and C.-W. Shu. Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: One dimensional case. J. Comput. Phys., 212:540–565, 2003. [44] J. Qiu and C.-W. Shu. A comparison of troubled cell indicators for Runge-Kutta discontinuous Galerkin methods using WENO limiters. SIAM J. Sci. Comput., 27:995–1013, 2005. [45] J. Qiu and C.-W. Shu. Hermite WENO schemes and their application as lim- iters for Runge-Kutta discontinuous Galerkin method ii: Two dimensional case. Comput. & Fluids, 34:642–663, 2005. [46] J. Qiu and C.-W. Shu. Runge-Kutta discontinuous Galerkin method using WENO limiters. SIAM J. Sci. Comput., 26:907–929, 2005. [47] W. Reed and T. Hill. Triangular mesh methods for the neutron transport equa- tion. Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [48] D. R. Scott and D. J. Stevenson. Magma ascent by porous flow. J. Geophys. Res., 91:9283–9296, 1986. 120 [49] C.-W. Shu. TVB uniformly high-order schemes for conservation laws. Math. Comp., 49:105–121, 1987. [50] C.-W. Shu. TVD time-discretizations. SIAM J. Sci. Statist. Comput., 9:1073– 1084, 1988. [51] C.-W. Shu. Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws. Technical Report ICASE Report No. 97-65, NASA Langley Research Center, 1997. [52] C.-W. Shu and S. Osher. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comput. Phys., 77:439–471, 1988. [53] N. H. Sleep. Segregation of magma from a mostly crystalline mush. Geol. Soc. Am. Bull, 85:1225–1232, 1974. [54] M. Spiegelman. Flow in deformable porous media, 1, simple analysis. J. Fluid. Mech., 247:39–63, 1993. [55] M. Spiegelman and P. B. Kelemen. Extreme chemical variability as a conse- quence of channelized melt transport. Geochem. Geophys. Geosyst., 4(7), 2001. [56] M. Spiegelman, P. B. Kelemen, and E. Aharonov. Causes and consequences of flow organization during melt transport: The reaction infiltration instability in compactible media. J. Geo. Res., 106:2061–2077, 2001. [57] E. Stolper. A phase diagram for mid-ocean ridge basalts: Preliminary results and implications for petrogenesis. Contrib. Mineral. Petrol., 74:13–27, 1980. [58] G. Szego. Orthogonal Polynomials. American Mathematical Society, 1981. [59] L. Trefethen. Spectral Methods in Matlab. SIAM Publishing, 2000.