An Implementation of the Finite Element Method for the Velocity-Current Magnetohydrodynamics Equations by K. Daniel Brauss A dissertation submitted to the Graduate Faculty of Auburn University in partial ful llment of the requirements for the Degree of Doctor of Philosophy Auburn, Alabama August 2, 2014 Keywords: magnetohydrodynamics, nite element method, numerical solution Copyright 2014 by K. Daniel Brauss Approved by A. J. Meir, Chair, Professor of Mathematics and Statistics Paul G. Schmidt, Professor of Mathematics and Statistics Yanzhao Cao, Professor of Mathematics and Statistics Dmitry Glotov, Professor of Mathematics and Statistics Abstract The aim of this research is to extend the numerical results in [78] of a velocity-current magnetohydrodynamics formulation proposed by A.J. Meir and Paul G. Schmidt in [74] for a stationary ow that models a conductive uid in a bounded domain. A parallel - nite element algorithm was successfully implemented on a high-performance computing, distributed-memory architecture at the Alabama Supercomputer Center (ASC) using the freely available, open-source academic and government libraries deal.ii, p4est and Trilinos. Extending the work of Elman, Silvester and Wathen [37] for the Navier-Stokes equations, a Schur complement preconditioner was developed for the current saddle-point problem to successfully utilize the iterative Krylov subspace solver GMRES (generalized minimal resid- ual method) and solve large linear systems of equations arising from mesh re nement. To simplify and lower operation costs in forming the preconditioner, spectral equivalence was established between the Schur complement and a mass matrix. The resulting C++ code was tested succesfully on problems from [78] with similar results. ii Acknowledgments I have to thank God for the opportunity to continue my studies here at Auburn University. I have to thank my wife Minvera Rosario Brauss Blanco for her unwavering support. Without my wife, I would not have made it this far. Thank you Rose. I thank Eliana Catherine Brauss, my daughter, for being here and lling our lives. I thank all of my family for their support with all that we do. I thank all of my teachers for the passion that they have given to me. I thank Dr. A. J. Meir for his much appreciated support, patience and guidance as my Ph.D. advisor. Dr. Meir is a positive in uence and a well-respected role model. I have learned a deeper independence from him that has helped me to grow. I thank Dr. Paul G. Schmidt for his support and his courses in functional analysis. Thanks to my other committee members Dr. Yanzhao Cao and Dr. Dmitry Glotov. Lastly, thanks to Dr. Narendra Govil and Dr. Geraldo DeSouza. iii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 Preconditioner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.4 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 The Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.1 Problem Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Weak Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3 Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3 Matrix System and Preconditioner . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.1 Matrix Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Finite Element Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3 Preconditioner . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 iv 4 Computational Experiments and Results . . . . . . . . . . . . . . . . . . . . . . 47 4.1 Introduction to the Computational Environment . . . . . . . . . . . . . . . . 47 4.1.1 Choice of the Libraries . . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.1.2 Working with the Libraries . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 Example 1 - A Problem with an Exact Solution . . . . . . . . . . . . . . . . 51 4.3 Example 2 - An Applied Problem . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.1 Preconditioning Results . . . . . . . . . . . . . . . . . . . . . . . . . 58 5 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.2 Torus Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.3 Code Speedup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 A Weak Formulation and Identities . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.1 Navier-Stokes Identities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.1.1 Pressure Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 A.1.2 Di usion Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 A.1.3 Advection Term . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 A.2 Ohms Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 B Newton?s Method Linearization . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 C Iterative Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 C.1 Conjugate Gradient (CG) Method . . . . . . . . . . . . . . . . . . . . . . . . 89 v C.2 Generalized Minimal Residual (GMRES) Method . . . . . . . . . . . . . . . 91 D Exact Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 vi List of Figures 2.1 Newton?s Method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.1 2-D Example of Q1 Basis Function for Pressure . . . . . . . . . . . . . . . . . . 30 3.2 2-D Example of Q2 Basis Function for Velocity . . . . . . . . . . . . . . . . . . 30 4.1 Schafer Turek Benchmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Example 1 for Cube Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3 Example 1 for Lshaped Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.4 Example 2 for Cube Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.5 Example 2 for Lshaped Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.6 Eigenvalues for System Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.7 Eigenvalues for Preconditioned Stokes-Ohms System . . . . . . . . . . . . . . . 60 4.8 Eigenvalues for Preconditioned Stokes-Ohms System with Inertial Term . . . . . 61 4.9 Eigenvalues for Preconditioned System Matrix . . . . . . . . . . . . . . . . . . . 61 vii 5.1 Example 1 for Lshaped Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Example 2 for Torus Domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.3 Example 2 for Full Torus Domain . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.4 CPU Time Versus Number of Processors . . . . . . . . . . . . . . . . . . . . . . 68 5.5 CPU Time Versus Number of Processors . . . . . . . . . . . . . . . . . . . . . . 70 viii List of Tables 4.1 Errors and Convergence Rates for the Unit Cube Domain . . . . . . . . . . . . 53 4.2 Errors and Convergence Rates for the Lshaped Domain . . . . . . . . . . . . . . 53 5.1 Errors and Convergence Rates for the Quarter Torus Domain . . . . . . . . . . 63 ix Chapter 1 Introduction 1.1 Overview Scienti c theories often have a quantitative, mathematical nature to them, such as a well- founded set of governing equations developed through experimental tests and results. The equations ususally depict a variable, quantity or quantities that are of scienti c interest and may be a direct measure of the success of a process or system (e.g., velocity in mass transfer) or an indirect measure of a system?s state (e.g., pressure or concentration pro les indicating compression or swelling). Models whose governing equations involve complex multiphysics and systems of partial di erential equations often arise in physics and engineering. Advances in scienti c computing make it possible to numerically approximate solutions to the quanti- ties of interest in such equations with increased accuracy and less simplifying assumptions. Numerically approximating and visualizing solutions to complex multiphysics models can then allow a scientist to investigate and analyze a process, assess its conditions, and address possible optimizations. With growth in computing, the area of computational science and applied mathematics has established itself as a the third pillar of scienti c investigation [31]. By way of visualization, scienti c computing can o er a deep understanding of phenomena being scienti cally investigated as well as the limits of our theory and understanding of natural phenomena [23]. In this research, our interests are in the predictions of scienti c quantities according to governing equations for electrically conducting uid ow. We are interested in the dynamics of electrically conductive uids in the presence of electromagnetic elds, placing our research 1 in the area of magnetohydrodynamics. We study viscous, incompressible, homogeneous, conductive uids and therefore the governing incompressible Navier-Stokes equations [69] @v@t u + (u r) u +rp = F + J B r v = 0: We consider the Lorentz force J B to act on the uid by way of electrical and magnetic elds governed by the Maxwell equations [51, 38] and Ohm?s Law r E = q 0 (Gauss?s Law for Electric Field); r B = 0 (Gauss?s Law for Magnetic Field); r E = @B@t (Faraday?s law); r B = 0 J + 0@E@t (Ampere-Maxwell law); and J = (E + u B) (Ohm?s Law) The quantities of interest are the current J, velocity u, pressure p, electric potential , and magnetic eld B. Magnetohydrodynamics (MHD) is the study of the dynamics of electrically conducting uids under electromagnetic elds. Our focus is on the interaction of electrically conducting uids with electromagnetic elds at the macroscopic level [29]. Applications modeled by such equations arise in liquid metal solidi cation processes [77, 83, 84, 41, 42, 70, 90, 49, 79, 32, 18] and silicon crystal growth processes used in the semiconductor and microelectronics 2 industry [50, 73, 56, 55] to hold, liquify and stir melts, control undesired convection, and lter out impurities. Other example applications are electromagnetic pumping, nuclear reactor cooling, plasmas in nuclear fusion and propulsion devices [66, 54, 76, 52, 62, 57, 3, 72, 71, 88]. We note that these studies address grand challenges seen in the National Academy of Engineering?s (NAE) 21st Century?s Grand Challenges for Engineering [82] and the 5 Grand Challenges from the Basic Energy Sciences Advisory Committee (BESAC) [39, 58] of the U.S. Department of Energy (DOE). One particular grand challenge from NAE is providing energy from thermonuclear fusion - the harnessing of the sun?s energy. The above research also addresses the main theme of the BESAC grand challenges - designing and controlling material processes. The velocity-current formulation of the MHD equations that we study was proposed by A.J. Meir and Paul G. Schmidt [74]. They have established existence and uniqueness of a solution to the equations as well as a convergent nite element iterative numerical approximation method [78]. The equations are a stationary (steady-state) form of the equations above and will be reviewed in chapter 2. The main goal of this research was to extend the work conducted by A.J. Meir and Paul G. Schmidt in the numerical solution of the velocity-current MHD equations on three- dimensional domains by implementing a nite element method using parallel computing and a high-performance distributed memory cluster. To provide some motivation for this e ort, we consider the ow of a uid contained in a simple domain such as the unit cube. If we partition the unit cube in each coordinate direction into 101 segments to create a 100-by-100-by-100 domain mesh, and we perform the nite element method on this mesh, then the resulting linear system of equations will have unknowns that reach well into the millions. We have interest to apply our method to more complicated and larger domains as well as extend our work to include further multiphysics. This would mean that the number of unknowns could reach into the billions. We note that we are solving a nonlinear system of the partial di erential equations using the nite element method and a picard iteration 3 scheme, and therefore solve this size of a system repeatedly to convergence. As might be expected, a large portion of computational time is spent solving the linear system and is common in scienti c computing [99, 94]. 1.2 Solver We therefore wish to set up and e ciently solve a linear system Ax = b. Two general choices for solving a linear system are a direct solver or an iterative solver. A direct solver can be viewed as a form of Gaussian elimination performed on the system to decompose it into lower L and upper U triangular factors, A = LU. Many iterative solvers can be viewed as projection methods [85, 95] where the iterative solution to the system is obtained by projecting the exact solution to the system onto a particular subspace using a technique such as least squares. As our linear system is large, our interest is in the operations and memory required for such solvers. We summarize from a serial point of view. The operation count for Gaussian elimination on a n n matrix A, excluding row inter- changes, is O(n3). Indeed, the pivot on the main diagonal element at the (1,1) position results in n multiplications across the rst row. The elimination (zeroing) of the element below the pivot in the rst column requires a multiplication on the rst row (n ops) and a row addition (n ops). Performing elimination for each of the n 1 rows below the pivot results in a total of n+ (n+n)(n 1) = n+ 2n(n 1) ops ( oating point operations). The total number of operations across each column can be written as a series n 1X i=0 (n i) + 2(n i)(n 1 i) = n 1X i=1 i+ 2(i)(i 1) = 2 n 1X i=1 i2 n 1X i=1 i = (n 1)(n)(2n 1)3 (n 1)(n)2 = O(n3): 4 To solve the system Ax = b, the same operations are applied to the system?s right-hand side b creating a column vector we call c and requiring n+ n 1X i=1 2(i 1) = n(n 1) n = O(n2) oating point operations. Upon completion, the system of equations has been transformed to Ux = c. To obtain the solution x, a backsolve is performed, starting from the lower row of U and working upward. The number of operations to solve the ith element of x is one multiplication and i 1 additions. The complexity for solving all the elements of x, the entire backsolve, is nX i=1 1 + nX i=1 (i 1) = nX i=1 i = n(n+ 1)2 = O(n2): The order of operations for the entire solution process is therefore O(n3). We note that bands in the system matrix are common for numerical approximation techniques to solutions of partial di erential equations, such as the nite element method, and this matrix structure is taken advantage of by direct solvers. For a banded matrix with b nonzero diagonals below the main diagonal, the complexity when solving with Gaussian elimiation is less and dependent on the band size O(nb2) [44]. To compare the number of operations between a direct and iterative solver, we take as an example a second order PDE discretized using the nite element method over an irregular three-dimensional grid to obtain the system of equations Ax = b. A typical bandwith for this problem is b = O(n2=3) [97]. The number of operations using the banded direct solver would therefore be O(nb2) = O(n7=3). For an iterative solver, the symmetry and positive de niteness of the system matrix makes the conjugate gradient (CG) method a natural choice. Per iteration the number of oating point operations is O(n), see appendix C. To 5 estimate the number of operations for the complete iterative solve, we use an error bound after l iterations of the method given in the following theorem [44]. Theorem 1.1 Suppose that A2Rn n is symmetric positive de nite and b2Rn. The error after l iterations of the CG algorithm in solving the system Ax = b can be bounded as follows: kx xlkA 2 p (A) 1p (A) + 1 !l kx x0kA; where (A) is the condition number of A. For discretized second order PDEs, it can be shown that (A) = O(n2=3) [97]. Setting the iterative solver error tolerance to , we estimate the required number of iterations l to achieve the tolerance by using the error estimate from the theorem above. Here we have assumed that n is large and used (A) = O(n2=3) to simplify. 1 1=p (A) 1 + 1=p (A) !l < 1 2p (A) !l