Iterative Solution of a Mixed FiniteElement Discretisation of anIncompressibleMagnetohydrodynamics ProblembyMichael WathenM.Sci., The University of Birmingham, 2012A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Computer Science)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2014c© Michael Wathen 2014AbstractThe aim of this thesis is to develop and numerically test a large scale pre-conditioned finite element implementation of an incompressible magnetohy-drodynamics (MHD) model. To accomplish this, a broad-scope code hasbeen generated using the finite element software package FEniCS and thelinear algebra software PETSc. The code is modular, extremely flexible,and allows for implementing and testing different discretisations and lin-ear algebra solvers with relatively modest effort. It can handle two- andthree-dimensional problems in excess of 20 million degrees of freedom.Incompressible MHD describes the interaction between an incompress-ible electrically charged fluid governed by the incompressible Navier-Stokesequations coupled with electromagnetic effects from Maxwell’s equations inmixed form. We introduce a model problem and a mixed finite elementdiscretisation based on using Taylor-Hood elements for the fluid variablesand on a mixed Ne´de´lec pair for the magnetic unknowns. We introducethree iteration strategies to handle the non-linearities present in the model,ranging from Picard iterations to completely decoupled schemes.Adapting and extending ideas introduced in [Dan Li, Numerical Solu-tion of the Time-Harmonic Maxwell Equations and Incompressible Magne-tohydrodynamics Problems, Ph.D. Dissertation, The University of BritishiiColumbia, 2010], we implement a preconditioning approach motivated bythe block structure of the underlying linear systems in conjunction withstate of the art preconditioners for the mixed Maxwell and Navier-Stokessubproblems. For the Picard iteration scheme we implement an inner-outerpreconditioner.The numerical results presented in this thesis demonstrate the efficientperformance of our preconditioned solution techniques and show good scal-ability with respect to the discretisation parameters.iiiPrefaceThis dissertation is original, unpublished, independent work by the author,M. Wathen.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 A model problem in incompressible magnetohydrodynamics . 11.2 Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . 61.2.1 Finite element methods for incompressible MHD prob-lems . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.2 Preconditioning incompressible MHD problems . . . . 91.3 Objectives and contributions . . . . . . . . . . . . . . . . . . 101.4 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11v2 Finite element discretisation . . . . . . . . . . . . . . . . . . . 132.1 Variational formulation . . . . . . . . . . . . . . . . . . . . . 132.2 Mixed finite element discretisation . . . . . . . . . . . . . . . 162.2.1 Matrix representation . . . . . . . . . . . . . . . . . . 182.3 Picard iteration (P) . . . . . . . . . . . . . . . . . . . . . . . 202.4 Decoupled iterations . . . . . . . . . . . . . . . . . . . . . . . 232.4.1 Magnetic Decoupling (MD) . . . . . . . . . . . . . . . 242.4.2 Complete Decoupling (CD) . . . . . . . . . . . . . . . 252.5 Inhomogeneous Dirichlet boundary conditions and initialguess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273 Preconditioning . . . . . . . . . . . . . . . . . . . . . . . . . . . 283.1 Preconditioning the incompressible Navier-Stokes equations . 293.1.1 Pressure Convection-Diffusion (PCD) . . . . . . . . . 303.1.2 Least-Squares Commutator (LSC) . . . . . . . . . . . 313.1.3 PCD versus LSC . . . . . . . . . . . . . . . . . . . . . 333.2 Preconditioning Maxwell’s equations . . . . . . . . . . . . . . 343.2.1 An ideal preconditioner . . . . . . . . . . . . . . . . . 343.2.2 A practical preconditioner . . . . . . . . . . . . . . . 353.3 Preconditioning the MHD equations . . . . . . . . . . . . . . 363.3.1 Picard iteration . . . . . . . . . . . . . . . . . . . . . 363.3.2 Magnetic Decoupling . . . . . . . . . . . . . . . . . . 383.3.3 Complete Decoupling . . . . . . . . . . . . . . . . . . 393.3.4 Summary . . . . . . . . . . . . . . . . . . . . . . . . . 40vi4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . 424.1 Software . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . 434.3 Numerical results: Navier-Stokes equations . . . . . . . . . . 454.3.1 Convergence results for a 2D smooth solution . . . . 464.3.2 Convergence results for a 3D smooth solution . . . . 474.3.3 Preconditioning with LSC and PCD . . . . . . . . . . 484.4 Numerical results: Maxwell’s equations . . . . . . . . . . . . 494.4.1 Convergence results for a 2D smooth solution . . . . 504.4.2 Convergence results for a 3D smooth solution . . . . 514.4.3 Preconditioning . . . . . . . . . . . . . . . . . . . . . 524.5 Numerical results: MHD problem . . . . . . . . . . . . . . . 544.5.1 Convergence results for a 2D smooth solution . . . . 544.5.2 Convergence results for a 3D smooth solution . . . . 564.5.3 Parameter tests . . . . . . . . . . . . . . . . . . . . . 584.6 Preconditioned MHD problem . . . . . . . . . . . . . . . . . 604.6.1 Picard iteration . . . . . . . . . . . . . . . . . . . . . 614.6.2 Magnetic Decoupling . . . . . . . . . . . . . . . . . . 634.6.3 Complete Decoupling . . . . . . . . . . . . . . . . . . 655 Conclusions and future work . . . . . . . . . . . . . . . . . . . 695.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 695.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76viiAppendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86A Curl operators and cross products . . . . . . . . . . . . . . . 87A.1 2D curl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87A.2 3D curl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88B Krylov subspace methods . . . . . . . . . . . . . . . . . . . . . 89viiiList of Tables3.1 Summary of coefficient matrices and corresponding precondi-tioners for each iteration scheme . . . . . . . . . . . . . . . . 404.1 Convergence of the velocity field for 2D Navier-Stokes -tol = 1e-10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.2 Convergence of the pressure variable for 2D Navier-Stokes -tol = 1e-10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 464.3 Convergence of the velocity field for 3D Navier-Stokes -tol = 1e105 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 474.4 Convergence of the pressure variable for 3D Navier-Stokes -tol = 1e-10 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 484.5 Iteration table for LSC preconditioner for 2D example forvarious values of ν and tol = 1e-5 . . . . . . . . . . . . . . . . 494.6 Iteration table for a PCD preconditioned 2D example for var-ious values of ν and tol = 1e-5 . . . . . . . . . . . . . . . . . 504.7 Convergence for 2D Maxwell smooth solution - magnetic field 504.8 Convergence for 2D Maxwell smooth solution - multipliervariable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 514.9 Convergence for 3D Maxwell smooth solution - magnetic field 52ix4.10 Convergence for 3D Maxwell smooth solution - multipliervariable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 524.11 Iteration count for Maxwell preconditioner for 2D example -direct application of preconditioner . . . . . . . . . . . . . . . 534.12 Iteration table for Maxwell preconditioner for 3D example -direct application of preconditioner . . . . . . . . . . . . . . . 534.13 Convergence for 2D MHD smooth solution (velocity field) -tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.14 Convergence for 2D MHD smooth solution (pressure variable)- tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.15 Convergence for 2D MHD smooth solution (magnetic field) -tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 554.16 Convergence for 2D MHD smooth solution (multiplier vari-able) - tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . 564.17 Convergence for 3D MHD smooth solution (velocity field) -tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.18 Convergence for 3D MHD smooth solution (pressure variable)- tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.19 Convergence for 3D MHD smooth solution (magnetic field) -tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 574.20 Convergence for 3D MDH smooth solution (multiplier vari-able) - tol = 1e-8 . . . . . . . . . . . . . . . . . . . . . . . . . 574.21 Number of non-linear iterations for various values of ν withtol = 1e-5, κ = 1 and νm = 10. . . . . . . . . . . . . . . . . . 59x4.22 Number of non-linear iterations for various values of κ withtol = 1e-5, ν = 1 and νm = 10. . . . . . . . . . . . . . . . . . 594.23 Number of non-linear and average number of preconditioningiterations for various values of ν with tol = 1e-5, κ = 1 andνm = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 624.24 Number of non-linear and average number of preconditioningiterations for various values of κ with tol = 1e-5, ν = 1 andνm = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 634.25 Number of non-linear iterations and average number of iter-ations to solve the Navier-Stokes and Maxwell’s subproblemfor the MD scheme with tol = 1e-4, κ = 1, ν = 1 and νm = 10. 634.26 Number of non-linear iterations and average number of iter-ations to solve the Navier-Stokes and Maxwell’s subproblemfor the MD scheme with tol = 1e-4, ν = 1 and νm = 10. . . . 644.27 Number of non-linear iterations and average number of iter-ations to solve the Navier-Stokes and Maxwell’s subproblemfor the MD scheme with tol = 1e-4, κ = 1, ν = 1 and νm = 10. 654.28 Number of non-linear iterations and average number of iter-ations to solve the Navier-Stokes and Maxwell’s subproblemfor the MD scheme with tol = 1e-4, κ = 1, ν = 1 and νm = 10in 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.29 Number of non-linear iterations and average number of iter-ations to solve the Stokes and Maxwell’s subproblem for theCD scheme with tol = 1e-4 κ = 1, and νm = 10. . . . . . . . . 66xi4.30 Number of non-linear iterations and average number of iter-ations to solve the Stokes and Maxwells subproblem for theCD scheme with tol = 1e-4, ν = 1 and νm = 10. . . . . . . . . 674.31 Number of non-linear iterations and average number of iter-ations to solve the Stokes and Maxwell’s subproblem for theCD scheme with tol = 1e-4, κ = 1, ν = 1 and νm = 10. . . . . 674.32 Number of non-linear iterations and average number of iter-ations to solve the Stokes and Maxwell’s subproblem for theCD scheme with tol = 1e-4, κ = 1, ν = 1 and νm = 10 in 3D. 68xiiList of Figures4.1 Level 3 grid on unit square domain . . . . . . . . . . . . . . . 44xiiiAcknowledgementsFirst, and foremost, I would like to thank my supervisors Chen Greif andDominik Scho¨tzau for their help and guidance. I greatly appreciate theirunderstanding and insights throughout my graduate study and writing ofthis thesis. I thank Uri Ascher for reading this thesis and providing me withexcellent feedback.I would like to thank UBC Computer Science for their financial supportover the last few year.Finally, I would like to thank my friends and family for there love andsupport throughout my studies.xivChapter 1IntroductionThe primary topic of this thesis is to develop and numerically test a largescale implementation of an incompressible magnetohydrodynamics model.In this introductory chapter, we will first present a description of the modelproblem studied. We then give a brief overview of finite element methodsand Krylov subspace methods for this problem. Finally, we outline theobjectives, contributions and structure of the thesis.1.1 A model problem in incompressiblemagnetohydrodynamicsThe area of incompressible magnetohydrodynamics (MHD) describes thebehaviour of electrically conductive incompressible fluids (liquid metals,plasma, salt water, etc) in an electromagnetic field [14, 23, 45]. MHD modelscouple electromagnetism and fluid dynamics. The coupling effects are dueto two fundamental physical properties. Firstly, through the movement ofthe conductive material that induces a magnetic field which then modifiesany existing electromagnetic field. Secondly, the magnetic and electric fieldsgenerate a mechanical force on the fluid known as the Lorentz force. The1Lorentz force accelerates the fluid particles in the direction normal to boththe electric and magnetic fields.Incompressible MHD has a number of important applications withintechnology and industry as well as Geophysical and Astrophysical appli-cations. Some such applications are: electromagnetic pumping, aluminiumelectrolysis, the Earth’s molten core and solar flares. For more applicationssee [45].In this thesis, we are principally interested in an incompressible MHDmodel. This means that the electrically conductive fluid is incompressible,i.e., the mass of the fluid is conserved, and the electric resistivity of thefluid cannot be ignored. The MHD model we consider consist of two cou-pled fundamental equations: the incompressible Navier-Stokes equationsand Maxwell’s equations. We will outline the derivation of a formulationof an incompressible MHD model for a homogeneous and isotropic mediumto ensure that all material parameters are constant; for full details see [8].The transient incompressible Navier-Stokes equations that govern in-compressible fluid flow are given by:ρf(∂u∂t+ (u · ∇)u)− µ∆u+∇p = f + fL in Ω× (0, T ), (1.1a)∇ · u = 0 in Ω× (0, T ). (1.1b)Here u and p are the velocity and pressure of the fluid, f denotes the bodyforces acting on the fluid and fL is the Lorentz force, which will be specifiedlater. The parameters µ > 0 and ρf > 0 denote the dynamic viscosity anddensity of the fluid, respectively. The spatial domain is given by Ω and the2end time is denoted by T . Mass conservation is given by (1.1b), see [18,Chapter 0] for the derivation of the incompressible Navier-Stokes equations.Maxwell’s equations that govern electromagnetic effects are given by:Faraday’s law: ∂b∂t+∇× e = 0, (1.2a)Coulomb’s law: ∇ · d = ρˆe, (1.2b)Ampe`re’s law: −∂d∂t+∇× h = j, (1.2c)Gauss’s law: ∇ · b = 0. (1.2d)The fields in (1.2) are given by: h the magnetic field, e the electric field, dthe electric displacement, b the magnetic induction and j the electric currentdensity. The parameter ρˆe is the electric charge density. In a homogeneousand isotropic medium, the following linear relations hold:d = δe, b = µh, (1.3)where the constant δ > 0 denotes the electric permittivity and the constantµ > 0 the magnetic permeability. Using (1.2) together with (1.3) yields theform of Maxwell’s equations considered in this thesis:Faraday’s law: ∂b∂t+∇× e = 0, (1.4a)Coulomb’s law: ∇ · e = ρe, (1.4b)Ampe`re’s law: −∂(δe)∂t+∇× ( 1µb) = j, (1.4c)Gauss’s law: ∇ · b = 0, (1.4d)3where ρe = ρˆeδ ; see [44, Chapter 1] for more details on Maxwell’s equations.The physical assumptions we consider to form a MHD model are thesame as in [8, Section 2.1]. More precisely we assume:• Non-relativistic motion: The characteristic fluid velocity is assumedto be orders of magnitude smaller than the speed of light.• Low-frequency approximation: Phenomena involving high frequenciesare omitted. That is, the term ∂(δe)∂t involving the displacement currentis neglected in Maxwell’s equations. Therefore, Ampe`re’s law (1.4c)simplifies to∇×(1µb)= j. (1.5)• Quasi-neutrality assumption: Positive and negative charges are equalin any given region. The convection current is omitted and Ohm’s lawnow readsj = θ(e+ u× b), (1.6)where positive parameter θ defines the electric conductivity of thefluid and u× b corresponds to the charge density induced by the fluidmotion. The electrical resistivity is given by 1/θ and causes dissipativeeffects in Maxwell’s equations but is not to be neglected in this model.• Non-magnetisation and non-polarisation: The assumptions of homo-geneity and isotropy also imply that the medium is non-magnetisableand non-polarisable.4The Lorentz force fL in (1.1a) can now be expressed as:fL =1ρf µ(∇× b)× b. (1.7)The electromotive term in (1.6) now enters Faraday’s law (1.4a) by combin-ing (1.6) and (1.5) into the following expression for e:e = 1θ(∇×(1µb)− u× b).Using the assumptions above, elimination of the electric field e and non-dimensionalisation, the systems in (1.1) and (1.4) are coupled into the fol-lowing set of partial differential equations:∂u∂t− ν ∆u+ (u · ∇)u+∇p− κ (∇× b)× b = f , (1.8a)∇ · u = 0, (1.8b)∂b∂t+ νm∇× (∇× b)− ∇× (u× b) = 0, (1.8c)∇ · b = 0, (1.8d)with suitable initial conditions and boundary conditions; see [8, Section 2].The unknowns u, p and b are the fluid velocity, the fluid pressure and themagnetic field, respectively.The solution to (1.8) depends on three non-dimensional parametersν = 1/Re, νm = 1/Rm and κ. The first parameter Re is the hydrodynamicReynolds number, which indicates the balance between the inertial forcesand the viscous forces. The parameter Rm is the magnetic Reynolds number,5which measures the effect by which the magnetic field induces flow motion.The final parameter, the coupling number, κ, represents the influence of theelectromagnetic field on the flow. It is sometimes defined in terms of theHartmann number denoted by Ha, asHa =√κReRm.To find typical physical values for these parameters, we refer to [8, 23, 51].We refer to ν as the viscosity for the rest of this thesis.In this thesis, we are interested in the steady-state ( ∂∂t = 0) versionof (1.8):−ν ∆u+ (u · ∇)u+∇p− κ (∇× b)× b = f , (1.9a)∇ · u = 0, (1.9b)κνm∇× (∇× b)− κ∇× (u× b) = 0, (1.9c)∇ · b = 0. (1.9d)Note we have multiplied (1.9c) by κ to enforce skew-symmetry of the cou-pling terms. We will consider both two- and three-dimensional solutions to(1.9). See Appendix A for curl definitions.1.2 Numerical solutionThe partial differential equation (PDE) system given in (1.9) requires anumerical approximation as in general an analytical solution is not possible.6There are two main components in computing a numerical solution of aPDE:1. Discretisation: take a continuous model and transfer it into a discretemodel;2. Solve: take the discretised model and solve for the unknowns.In this thesis, we use mixed finite element methods for discretising an MHDmodel problem, and solve it using preconditioned Krylov subspace methods.In the sequel, we briefly describe these components in the context of theapproach we take.1.2.1 Finite element methods for incompressible MHDproblemsThere are several finite element methods for discretising MHD problemsas in (1.9). A common approach in the literature is to approximate themagnetic field using standard nodal H1-conforming elements [8, 22, 28, 50].Such a formulation enables one to use the following vector calculus identity−∆b = ∇× (∇× b)−∇(∇ · b). (1.10)Since b is divergences-free it is then possible to apply an augmentationtechnique to replace the curl-curl operator with a vector Laplacian. Thisthen reduces one of the principal computational difficulties, namely the largenull-space of the curl-curl operator in (1.9c). However, one of the mainproblems using H1-conforming elements for the magnetic field is that for7non-convex domains (such as the 2D L-shaped domain with a reentrantcorner, or the 3D Fichera corner domain with reentrant edges and corners)the magnetic field will converge to a solution that is not correct around thesingular point, see [12, 13]. Therefore, we consider a mixed discretisationthat captures singular solutions correctly. One such family of elements areH(curl) conforming Ne´de´lec elements [47].To enable the use of Ne´de´lec elements for the magnetic field we usethe mixed formulation in [25, 55]. This leads to the following governingequations in a domain Ω:−ν ∆u+ (u · ∇)u+∇p− κ (∇× b)× b = f in Ω, (1.11a)∇ · u = 0 in Ω, (1.11b)κνm∇× (∇× b) +∇r − κ∇× (u× b) = g in Ω, (1.11c)∇ · b = 0 in Ω, (1.11d)where we have introduced the Lagrange multiplier, r, in the form of ∇rin (1.11c). Again, u and p are the velocity and pressure of the fluids andb is the magnetic field. The introduction of the Lagrange multiplier, r,corresponds to the divergence-free constraint (1.11d) of the magnetic field.With the addition of the Lagrange multiplier, r, we may also introduce ageneric forcing term g associated with the Maxwell part of (1.11).The numerical tests that we will consider will have inhomogeneousDirichlet boundary conditions for the fluid and magnetic fields and homo-8geneous Dirichlet boundary condition for the multiplier, r, of the form:u = uD on ∂Ω, (1.12a)n× b = n× bD on ∂Ω, (1.12b)r = 0 on ∂Ω, (1.12c)where uD and bD are given functions, and n is the unit outward normal tothe boundary ∂Ω. Notice, by taking the divergence of equation (1.11c) weobtain Poisson’s equation (as ∇ · ∇ × b = 0):∆r = ∇ · g in Ω, r = 0 on ∂Ω.In many physical applications g is divergence-free, which implies that themultiplier, r, is zero. In general, the main purpose of the magnetic multiplieris to provide stability, see [17].1.2.2 Preconditioning incompressible MHD problemsIncompressible MHD problems have been extensively studied in the contextof various discretisation and formulations. However, the development ofpreconditioned iterative solutions to MHD problems is limited. We referthe reader to the Appendix for a review of Krylov subspace solvers forlinear systems. In this thesis our focus will be on the implementation ofpreconditioners which are tailored to MHD model problem.In the literature there have not been too many approaches to precodition-ing the MHD equations given in either the non-multiplier form (1.9) or with9the multiplier form (1.11). In the very recent work [50], an operator-basedpreconditioner for the non-multiplier MHD equations (1.9) has been pro-posed in the context of H1-conforming elements for the magnetic field. Toform their preconditioner the authors use the identity (1.10) and a discretecommutator idea to form approximations to the Schur complements. Thisis based on an approach for a preconditioner to the Navier-Stokes systemas in [18, Chapter 8]. The preconditioners we employ are based on similarNavier-Stokes preconditioners but rely on the Maxwell preconditioner in [27]for H(curl) elements.1.3 Objectives and contributionsThe aim of this thesis is to develop and test fully scalable iterative solu-tion methods for the incompressible MHD model (1.11), (1.12) using nat-ural Taylor-Hood elements [59] for the fluid variables and Ne´de´lec mixedelement [47] pair for the magnetic variables. Our numerical results showgood scalability with respect to the mesh size. We provide several tests tostudy the performance of the preconditioners with respect to the relevantnon-dimensional parameters. We also present two- and three-dimensionalresults.To enable large scale preconditioned tests of the MHD model we use thefinite element software FEniCS [40] together with the linear algebra softwarefrom PETSc [9, 10]. Using these two principal software packages, experimentswere run in excess of 20 million degrees of freedom. As well, it provides anexample of how significant physical problems described by partial differential10equations can be solved by combining state of the art numerical softwarepackages. The aim is to release the code for public use.As stated above, little has been done in terms of preconditioned iterationsfor the MHD model other than the works of [50] for the non-multiplier form(1.9). Our approach is based on H(curl) elements for the magnetic field [55],and motivated by the preliminary results of [35] in the context of exactlydivergence-free elements for the velocity field. It combines preconditionersfor the incompressible Navier-Stokes and Maxwell’s equations in [18, 27, 36].The availability of our large-scale solvers and code will hopefully allowmore development and research into such MHD models.1.4 OutlineThis thesis is made up of five chapters and is structured as follows. InChapter 2, we introduce a mixed finite element approximation to the MHDsystem (1.11)-(1.12). The mixed approximation is based on a standard nodalTaylor-Hood finite element approximation for the velocity field and the pres-sure, together with a mixed Ne´de´lec element approximation for the magneticfield and the multiplier. Using this approximation, we introduce three pos-sible non-linear iteration schemes.In Chapter 3, we present an overview of the preconditioning ap-proaches for the individual subproblems separately, namely the incompress-ible Navier-Stokes and Maxwell’s equations. We then apply these precondi-tioning techniques to propose preconditoning strategies for the linearisationswhich arise from the three non-linear iteration schemes from Chapter 2. In11particular, for the coupled linearised Picard scheme we propose an inner-outer preconditioning approach.In Chapter 4, to perform numerical experiments in both two and threespatial dimensions we use the following two main software packages FEniCS[40] and PETSc [9, 10]. We show convergence results for the linearisedMHD system along with the incompressible Navier-Stokes and Maxwellsubproblems in isolation. Along with the convergence results we numeri-cally test the preconditioning approaches for the three non-linear iterationschemes, providing heuristic tests with respect to the dimensionless param-eters (ν, νm and κ; see (1.11)) and mesh size. These tests examine therobustness of both the preconditioners and the iteration schemes.Chapter 5 provides conclusions and outlines possible extensions for fu-ture work.12Chapter 2Finite element discretisationIn this chapter we introduce a mixed finite element discretisation for thesteady-state incompressible MHD problem (1.11), (1.12) that models elec-trically conductive fluids under the influence of a magnetic field. Followingthe setting in [55], we use curl-conforming elements for the magnetic fieldand conforming continuous elements for the velocity field. The resultingdiscretisation is verified though a series of numerical experiments which ap-pear later in Chapter 4. For simplicity, we initially only discuss in detailhomogeneous Dirichlet boundary conditions, that isu = 0 and n× b = 0. (2.1)Inhomogeneous conditions as in (1.12) are discussed in Section 2.5.2.1 Variational formulationSuppose that the domain Ω is a Lipschitz domain of Rd for d = 2, 3. Toexpress the problem (1.11), (1.12) in weak form we follow [55] and denotethe L2-inner product on L2(Ω)d by (·, ·)Ω, for d = 2, 3. We introduce the13standard Sobolev spacesV = H10 (Ω)d ={u ∈ H1(Ω)d : u = 0 on ∂Ω},Q = L20(Ω) = {p ∈ L2(Ω) : (p , 1)Ω = 0},C = H0(curl; Ω) ={b ∈ L2(Ω)d : ∇× b ∈ L2(Ω)d¯, n× b = 0 on ∂Ω},S = H10 (Ω) = {r ∈ H1(Ω) : r = 0 on ∂Ω}.(2.2)where d¯ = 2d− 3 is used to cover the 2D and 3D cases. We write ‖ · ‖L2(Ω),‖ · ‖H1(Ω) and ‖ · ‖H(curl;Ω) for the associated natural norms. More precisely,for vector fields u, b and a scalar function r the norms are defined as follows:‖u‖L2(Ω) =(∫Ωu · u dx)12,‖u‖H1(Ω) =(‖u‖2L2(Ω) + ‖∇u‖2L2(Ω))12 ,‖b‖H(curl,Ω) =(‖b‖2L2(Ω) + ‖∇ × b‖2L2(Ω))12 ,‖r‖L2(Ω) =(∫Ωr2 dx)12,‖r‖H1(Ω) =(‖r‖2L2(Ω) + ‖∇r‖2L2(Ω))12 ,where ‖∇u‖2L2(Ω) is given by:‖∇u‖2L2(Ω) =∫Ωd∑i,j=1(∇u)ij(∇u)ij dx12.14The weak formulation of the incompressible MHD system (1.11), (1.12) con-sists in finding (u, p, b, r) ∈ V ×Q×C × S such thatA(u,v) +O(u;u,v) + C(b;v, b) +B(v, p) = (f ,v)Ω, (2.3a)B(u, q) = 0, (2.3b)M(b, c)− C(b;u, c) +D(c, r) = (g, c)Ω, (2.3c)D(b, s) = 0, (2.3d)for all (v, q, c, s) ∈ V × Q × C × S. The individual variational forms aregiven byA(u,v) =∫Ων∇u : ∇v dx,O(w;u,v) =∫Ω(w · ∇)u · v dx,B(u, q) = −∫Ω(∇ · u) q dx,M(b, c) =∫Ωκνm(∇× b) · (∇× c) dx,D(b, s) =∫Ωb · ∇s dx,C(d;v, b) =∫Ωκ (v × d) · (∇× b) dx,(2.4)where ∇u : ∇v is defined as∇u : ∇v =d∑i,j=1(∇u)ij(∇v)ij .In [55] it has been shown that this formulation of the problem is discretelyenergy-stable and has a unique solution for small data (i.e. for small ν, νm,15κ and forcing terms f and g with small L2-norms).2.2 Mixed finite element discretisationConsider the domain Ω to be divided up into a regular and quasi-uniformmesh Th = {K} consisting of triangles (d = 2) or tetrahedra (d = 3) withmesh size h. Based on the function spaces defined in (2.2), our finite elementapproximation will be sought in the finite dimensional spaces given by:V h = {u ∈ H1(Ω) : u|K ∈ Pk(K)d, K ∈ Th },Qh = { p ∈ L2(Ω) ∩H1(Ω) : p|K ∈ Pk−1(K), K ∈ Th },Ch = { b ∈ H0(curl; Ω) : b|K ∈ Pk−1(K)d ⊕Rk(K), K ∈ Th },Sh = { r ∈ H10 (Ω) : r|K ∈ Pk(K), K ∈ Th },(2.5)for k ≥ 2. We define Pk(K) as the space of polynomials of total degree atmost k on K and Rk(K) as the space of homogeneous vector polynomials oftotal degree k on K that are orthogonal to the position vector x. Here wenote that we are using Pk/Pk−1 Taylor-Hood elements for the fluid unknowns(u, p) [59]. For the magnetic variables (b, r) we use the curl-conformingNe´de´lec element pair of the first kind [47]. These choices of finite elementsspaces V h, Ch, Qh and Sh imply that we have conforming subspaces toour Sobolev spaces V , C, Q and S, respectively. Then the finite elementsolution to (2.3) consists in finding (uh, ph, bh, rh) ∈ V h × Qh × Ch × Sh16such thatA(uh,v) + O˜(uh;uh,v) + C(bh;v, bh) +B(v, ph) = (f ,v), (2.6a)B(uh, q) = 0, (2.6b)M(bh, c)− C(bh;uh, c) +D(c, rh) = (g, c), (2.6c)D(bh, s) = 0, (2.6d)for all (v, q, c, s) ∈ V h ×Qh ×Ch × Sh.The forms A,M,B,D and C stay the same as on the continuous level.However, for the convection term O˜(·; ·, ·) we modify the form O(w;u,v) ina standard fashion to ensure the energy-stability propertyO˜(w;u,u) = 0, ∀w,u ∈ V h. (2.7)To do so we integrate by parts the convection form O(w;u,u) to obtain∫Ω(w · ∇)u · u dx =− 12∫Ω∇ ·wu · u dx+ 12∫∂Ωw · n|u|2 ds,recalling that n is the unit outward normal on ∂Ω. Therefore, we choosethe modified convection form O˜(w;u,v) asO˜(w;u,v) =∫Ω(w · ∇)u · v dx+ 12∫Ω∇ ·wu · v dx− 12∫∂Ωw · nu · v ds.By construction, property (2.7) is now satisfied. Note also that for homoge-neous boundary conditions as assumed in (2.1), the boundary integral term17in O˜ can be omitted.Again in [55] it has been shown that this variational formulation of aMHD problem is discretely energy-stable and has a unique solution forsmall data. Also, optimal order error estimates in the mesh size h havebeen derived for small data using the stability property (2.7). Namely, forsufficiently smooth solutions, we have the error bound‖u− uh‖H1(Ω) + ‖b− bh‖H(curl;Ω) + ‖p− ph‖L2(Ω) + ‖r − rh‖H1(Ω) ≤ Chk,for a constant C > 0 independent of the mesh size. In addition, the L2-norm error for the velocity field is of order O(hk+1) (as V h consists of a fullpolynomial space on each element). However, we cannot expect L2-normerrors of order O(hk+1) for the magnetic field (as Ch does not consist of afull polynomial space on each element).2.2.1 Matrix representationThe variational formulation (2.6) now can be converted into a matrix repre-sentation. To do this, we introduce the basis function for the finite elementspaces in (2.5):V h = span〈ψj〉nuj=1, Qh = span〈αi〉mui=1, (2.8)Ch = span〈φj〉nbj=1, Sh = span〈βi〉mbi=1. (2.9)The aim now is to find the coefficient vectors u = (u1, . . . , unu) ∈ Rnu , p =(p1, . . . , pmu) ∈ Rmu , b = (b1, . . . , bnb) ∈ Rnb , and r = (r1, . . . , rmb) ∈ Rmb18of the finite element functions (uh, ph, bh, rh) in terms of the chosen bases.As usual, this is done by writing the bilinear forms in (2.6) in terms of thefollowing stiffness matrices and load vectors:Ai,j = A(ψj ,ψi), 1 ≤ i, j ≤ nu,Bi,j = B(ψj , αi), 1 ≤ i ≤ mu, 1 ≤ j ≤ nu,Di,j = D(φj , βi), 1 ≤ i ≤ mb, 1 ≤ j ≤ nb,Mi,j = M(φj ,φi), 1 ≤ i, j ≤ nb,fi = (f ,ψi)Ω, 1 ≤ i ≤ nu,gi = (g,φi)Ω, 1 ≤ i ≤ nb.For the two non-linear forms, O˜ and C, we define the corresponding stiffnessmatrices with respect to given finite element functions wh ∈ V h and dh ∈Ch in the first argument and their associated coefficient vectors w and d asO(w)i,j = O˜(wh;ψj ,ψi), 1 ≤ i, j ≤ nu,C(d)i,j = C(dh;ψj ,φi), 1 ≤ i ≤ nb, 1 ≤ j ≤ nu.Thus, the numerical solution to (1.11) consists in solving the non-linear19systemA+O(u) BT CT (b) 0B 0 0 0−C(b) 0 M DT0 0 D 0upbr=f0g0, (2.10)where the vectors u ∈ Rnu , p ∈ Rmu , b ∈ Rnb , and r ∈ Rmb are the unknowncoefficients of the finite element functions.2.3 Picard iteration (P)The discrete system (2.10) is non-linear, and therefore appling a non-linearsolver to this problem is necessary. A common choice to deal with the non-linearity within the incompressible Navier-Stokes equations in isolation isto perform Oseen or Picard iterations [18]. This involves linearising aroundthe current velocity and solving for updates.For simplicity we only consider the linearly convergent Picard itera-tions. Since we have modified the convection form to be discretely energy-stable as in (2.7), an advantage of this approach is that the discreteconvection-diffusion operator is real positive. Thus, with small data thefixed-point/Picard iteration is a contraction. A more efficient non-linearsolver is Newton’s method which converges quadratically near the solution.However, applying Newton’s method is more involved, as it requires con-struction and solving linear systems associated with a Jacobian as well asfinding an initial guess sufficiently close to the solution. We leave the im-20plementation of Newton’s method or other non-linear solvers as an area ofpossible future work (Section 5.2).We adapt the fixed-point/Picard iterations to an approach for the fullMHD system, where we linearise around the current velocity and mag-netic fields. Given a current iterate (uh, ph, bh, rh) we solve for updates(δuh, δph, δbh, δrh) and introduce the next iterate by setting:uh→ uh + δuh, ph → ph + δph,bh → bh + δbh, rh → rh + δrh.In variational form, the updates (δuh, δph, δbh, δrh) ∈ V h ×Qh ×Ch × Share found by solving the Picard system (P):A(δuh,v) + O˜(u; δuh,v) + C(bh;v, δuh) +B(v, δph) = Ru(uh, bh, ph;v),B(δuh, q) = Rp(uh; q),M(δbh, c) +D(c, δrh)− C(bh; δuh,v) = Rb(uh, bh, rh; c),D(δbh, s) = Rr(bh; s),for all (v, q, c, s) ∈ V h ×Qh ×Ch × Sh. Note that this system is linearisedaround (uh, bh). The right-hand side linear forms correspond to the residual21at the current iteration (uh, ph, bh, rh) defined by:Ru(uh, bh, ph;v) = (f ,v)Ω −A(uh,v)− O˜(uh;uh,v)− C(bh;v, bh)−B(v, ph),Rp(uh; q) = −B(uh, q),Rb(uh, bh, rh; c) = (g, c)Ω −M(bh, c) + C(bh;uh, c)−D(c, rh),Rr(bh; s) = −D(bh, s),for all (v, q, c, s) ∈ V h ×Qh ×Ch × Sh.In [55] it is shown that for small data the Picard iteration (P) will con-verge to the exact solution for any initial guess.To formulate the variational form of the Picard iteration (P) in matrixform, let (u, p, b, r) be the coefficient vectors associated with (uh, ph, bh, rh)and (δu, δp, δb, δr) be the coefficient vectors of (δuh, δph, δbh, δrh). Then itcan readily seen that the Picard iteration (P) amounts to solving the matrixsystemA+O(u) BT C(b)T 0B 0 0 0−C(b) 0 M DT0 0 D 0δuδpδbδr=rurprbrr, (2.11)22withru = f −Au−O(u)u− C(b)T b−BT p,rp = −Bu,rb = g −Mu+ C(b)b−DT r,rr = −Db.(2.12)At each non-linear iteration, the right hand side vectors and matrices O(u)and C(b) in (2.12) and (2.11) respectively must be assembled with the so-lution coefficient vectors (u, p, b, r) of the current iterate. Here, the matrixA is symmetric positive definite (SPD), O(u) is non-symmetric and −C(b),C(b)T appear in a skew symmetric fashion. We also note that M is sym-metric positive semidefinite (SPSD) with nullity mb corresponding to thedimension of the scalar space Sh giving rise to the discrete gradients, see[27].2.4 Decoupled iterationsThe full MHD system (1.11), (1.12) is a coupled system consisting of theincompressible Navier-Stokes and Maxwell’s equations, coupled through thenon-linear skew symmetric coupling term C(b). In addition, the convectionterm O(u) is non-linear as well. These two terms make the numerical so-lution challenging. Therefore, if one or both of these terms is small thenit may be possible to iterate explicitly. In particular if the coupling term,C(b), is small then we may completely decouple the system into an incom-pressible Navier-Stokes problem and a Maxwell problem. The two resultingdecoupling schemes are what we call Magnetic and Complete Decoupling23and are both described below. Note that unlike the Picard iteration, thereis no small data guarantee that such iterations based on these decouplingschemes will converge; although we see convergence for reasonable values ofthe non-dimensional parameters.2.4.1 Magnetic Decoupling (MD)Consider first the situation where there is weak coupling within the sys-tem, that is when C(b) is small. Then it may be possible to drop theseterms to completely decouple the system into the two subproblems, theincompressibleNavier-Stokes and Maxwell’s equations. We will call this ap-proach Magnetic Decoupling (MD). Then the system(2.11) simplifies toA+O(u) BT 0 0B 0 0 00 0 M DT0 0 D 0δuδpδbδr=rurprbrr, (2.13)withru = f −Au−O(u)u− C(b)T b−BT p,rp = −Bu,rb = g −Mu+ C(b)b−DT r,rr = −Db.24We iterate in the same fashion as the Picard iteration with the simplermatrix (2.13). From (2.13) we can see that the system is now completelydecoupled. This enable us to solve each individual subproblem separatelyand possibly in parallel.2.4.2 Complete Decoupling (CD)For the second decoupling scheme, we again consider there to be weak cou-pling of the system but we also consider that the fluid equations are diffusion-dominated and hence can exclude the convection terms. This amounts tothe systemA BT 0 0B 0 0 00 0 M DT0 0 D 0δuδpδbδr=rurprbrr, (2.14)withru = f −Au−O(u)u− C(b)T b−BT p,rp = −Bu,rb = g −Mu+ C(b)b−DT r,rr = −Db.Again, we perform iterations in the same fashion as the Picard iteration.This is the simplest technique as it removes all non-linear terms in theiteration matrix and hence leaves the linear Stokes problem in the upper25(1, 1) block matrix.2.5 Inhomogeneous Dirichlet boundaryconditions and initial guessWhen considering inhomogeneous Dirichlet boundary conditions as in (1.12),we still solve (2.11), (2.13) and (2.14) for the solution updates with homoge-neous Dirichlet boundary conditions. Therefore, in this approach we mustincorporate the inhomogeneous Dirichlet boundary conditions only withinthe initial guess.To form a suitable initial guess, we solve the decoupled Stokes problemwith the inhomogeneous boundary condition (1.12a):A BTB 0up=f0, (2.15)and then the non-symmetric Maxwell problem the with inhomogeneousboundary conditions (1.12b), (1.12c):M − C(u) DTD 0br=g0. (2.16)Here the term C(u) corresponds the coupling term using u (the initial guessfor the velocity field). We expect the inclusion of the coupling term toincrease the accuracy of the initial guess because additional information ofthe problem is used.26The inhomogeneous Dirichlet boundary conditions are incorporated ina standard fashion by suitably modifying the matrix system. The outcomeof this procedure is that the boundary data interpolation is only performedfor the initial guess. Hence, the iterative solves for the initial guess must berun to a sufficient accuracy to ensure the accuracy of the discrete boundaryconditions.2.6 SummaryIn this chapter we reviewed a mixed finite element approximation to thefull MHD system given in (1.11) and (1.12). We followed the mixed ap-proach outlined in [55] and expressed the MHD system in the matrix form(2.11). Using the Picard iteration (2.11) we introduced two possible decou-pling schemes ((MD) and (CD)) which may be simpler to solve for. Theperformance of the resulting three non-linear iteration schemes depends onthe values of the parameters κ, ν and νm. The next chapter will discusspossible preconditioning approaches for these systems.In the sequel, we shall omit the dependence of O(u) and C(b) on u and b,respectively, and simply write O and C.27Chapter 3PreconditioningThe linear system (2.11) is typically sparse and of large dimensions, henceto efficiently solve for it we use a preconditioned iterative approach as pro-posed in [35]. We start by reviewing some preconditioning strategies for theincompressible Navier-Stokes and Maxwell subproblems in isolation. Fromthese techniques we will then introduce and numerically test preconditionersfor the full MHD system in linearised form as in (2.11).Let us introduce a number of matrices in addition to the ones alreadydefined in Chapter 2: the velocity mass matrix Q, the pressure mass ma-trix Qp, the pressure convection diffusion matrix Fp, the pressure Laplacianmatrix Ap, the scalar Laplacian matrix L and the mass matrix X byQi,j =∫Ωψj ·ψi dx, 1 ≤ i, j ≤ nu, (3.1a)(Qp)i,j =∫Ωαjαi dx, 1 ≤ i, j ≤ mu, (3.1b)(Fp)i,j = ν∫Ω∇αj · ∇αi + (w · ∇αj)αi dx, 1 ≤ i, j ≤ mu, (3.1c)(Ap)i,j =∫Ω∇αj · ∇αi dx, 1 ≤ i, j ≤ mu, (3.1d)Li,j =∫Ω∇βj · ∇βi dx, 1 ≤ i, j ≤ mb (3.1e)Xi,j =∫Ωφj · φi dx, 1 ≤ i, j ≤ nb, (3.1f)28where w is the finite element velocity from the current iteration. The matri-ces Fp and Ap are well defined as we use continuous elements for the pressurefinite element space Qh. The matrices incorporate homogeneous Neumannboundary conditions.3.1 Preconditioning the incompressibleNavier-Stokes equationsConsider the steady state incompressible Navier-Stokes equations in isola-tion. LetKNS =F BTB 0, (3.2)be the discretised and linearised Navier-Stokes (or Oseen) subproblem whereF = A + O has been defined in (2.4). Due to the convection term, O, thissystem is non-symmetric and we will use GMRES to iteratively solve thissubproblem; [54]. A common approach for solving a saddle point system, asin (3.2), is to use a block diagonal or block triangular preconditioner of theformMtNS =F BT0 −Sor MdNS =F 00 S, (3.3)where S ≈ BF−1BT is an approximation to the Schur complement. IfS = BF−1BT is precisely the Schur complement, then it has been provedin [46] that the preconditioned matrix has exactly 2 eigenvalues when usingthe block triangular preconditioner (i.e., M−1tNSKNS) or 3 distinct eigenvaluesin the block diagonal case (i.e., M−1dNSKNS). Both are diagonalisable and29hence GMRES will converge in exactly 2 or 3 iterations in the absence ofround-off errors.In practice it is often too expensive to form and solve for the exactSchur complement BF−1BT , hence, a good approximation is needed. Twowell-known preconditioners for the incompressible Navier-Stokes equationsare the Least Squares Commutator (LSC) and the Pressure Convection-Diffusion (PCD) preconditioners. A description of both can be found in[18] and we will just outline the procedure how these can be applied on thediscrete level. For the Navier-Stokes preconditioner we will consider blocktriangular preconditioners of the same form as MtNS.3.1.1 Pressure Convection-Diffusion (PCD)In [18, Chap. 8] the discrete commutator of the convection-diffusion operatorassociated with the gradient operation is introduced and given byh = (Q−1F )(Q−1BT )− (Q−1BT )(Q−1p Fp). (3.4)where the matrices Q, Qp and Fp are defined in (3.1). Assuming that thecommutator is small then pre- and post-multiplying (3.4) by BF−1Q andF−1p Qp, respectively, let us separate out the Schur complement to giveBF−1BT ≈ BQ−1BTF−1p Qp. (3.5)Our discretisation is inf-sup stable which implies that there is spectralequivalence between BQ−1BT and the pressure Laplacian matrix, Ap; see30[18, Section 5.5.1]. Hence, the Schur complement can be approximated by:SPCD = ApF−1p Qp.Applying the PCD preconditioner (i.e., taking S = SPCD) to the linearisedNavier-Stokes system involves solving the systemF BT0 −ApF−1p Qpxy=abat each Krylov iteration. This can be solved efficiently by splitting it intothe following two steps1. Solve for y: y = −Q−1p FpA−1p b;2. Solve for x: x = F−1(a−BT y).This means that we have one pressure Poisson solve (A−1p ), one mass ma-trix solve (Q−1p ), one convection-diffusion solve (F−1) and multiplicationswith Fp and BT at each Krylov iteration. It is possible to solve these it-eratively using multigrid and/or iterative methods. However, to test thepreconditioner we will use direct solver in this thesis.3.1.2 Least-Squares Commutator (LSC)One disadvantage with using the PCD preconditioner is that it requires theconstruction of the matrices Ap, Qp and Fp in (3.1). A second approachto approximate the Schur complement is the LSC preconditioner [18, Chap.8] which primarily uses the available matrix coefficients in (3.2) and the31construction of Q to form the preconditioner (without explicitly formingAp, Qp and Fp).As for the derivation of the PCD preconditioner we start off with thediscrete commutator of the convection-diffusion operatorh = (Q−1F )(Q−1BT )− (Q−1BT )(Q−1p Fp).Suppose that the Q-norm is defined by ‖v‖Q = (Qv, v)1/2. Then this timewe minimise h over the jth column of Fp (that is [Fp]j) in the Q-norm totry to find an approximation for Fp. As shown in [18], the minimisation isgiven bymin ‖[Q−1FQ−1BT ]j −Q−1BTQ−1p [Fp]j)‖Q.Solving this optimisation problem, as shown in [18], is done by solving thenormal equationsQ−1p BQ−1BTQ−1p Fp = Q−1p BQ−1FQ−1BT .This yields the an approximation to Fp asFp ≈ Qp(BQ−1BT )−1(BQ−1FQ−1BT ).By substituting this into expression (3.5) we obtain the LSC approximationto the Schur complement:S ≈ BF−1BT ≈ SLSC = (BQ−1BT )(BQ−1FQ−1BT )−1(BQ−1BT ).32Therefore, applying the LSC preconditioner to the Oseen system KNS in(3.2) involves solving for the matrixF BT0 −SLSCxy=abat each Krylov iteration. Again, this can be efficiently split up into thefollowing two steps:1. Solve for y: y = −(BQ−1BT )−1(BQ−1FQ−1BT )(BQ−1BT )−1b;2. Solve for x: x = F−1(a−BT y).Hence, we have two pressure Poisson solves ((BQ−1BT )−1) and oneconvection-diffusion solve (F−1) at each Krylov iteration as well as matrixmultiplications. In practice, we take the diagonal or lumped diagonal of Qto form BQ−1BT . These solves, as with the PCD preconditioner, will bedone directly in this thesis.3.1.3 PCD versus LSCThe main advantage of solving the commutator using the least-squares ap-proach is that the matrices that define the preconditioner are available fromthe original system KNS in (3.2) and the construction of Q . However, forPCD we require the construction of the matrices Ap, Qp and Fp. Therefore,LSC is slightly more computationally efficient to form. On the other hand,to apply the LSC preconditioner we require two pressure Poisson solves,whereas PCD only requires one.33We will consider experiments with both preconditioners for the incom-pressible Navier-Stokes problem in isolation to determine which seems moreeffective. This preconditioner will then be applied within the solver for thelinearised MHD system.3.2 Preconditioning Maxwell’s equationsNext, consider the Maxwell subproblemKMX =M DTD 0, (3.6)appearing in a linearised MHD system (2.11). As we did for the Navier-Stokes subproblem in Section 3.1, we apply a block preconditioning strategyfor KMX in (3.6). Notice that, KMX in (3.6) is symmetric and hence we willfocus on SPD block diagonal preconditioners.3.2.1 An ideal preconditionerThe (1, 1) block of KMX is the curl-curl operator, hence, the matrix Mis singular, where the null space is of dimension mb and consists of thediscrete gradients. Therefore the usual Schur complement does not exist asthe matrix M cannot be inverted. To overcome this difficulty, we employthe approach in [24, 26] based on augmentation and use preconditioners ofthe form:M +DTW−1D 00 W,34where W is a symmetric positive definite matrix.It has been shown in [27] that an appropriate choice of W is the scalarLaplacian, L, defined in (3.1). This leads to the ideal preconditioner:MiMX =M +DTL−1D 00 L. (3.7)Applying (3.7) as the preconditioner yields exactly two eigenvalues, 1 and−1 of algebraic multiplicities of nb and mb, respectively. Therefore usingthis matrix as a preconditioner means that MINRES will converge in twoiterations, in the absence of roundoff errors [49]. However, forming thematrix M+DTL−1D is costly, hence, MiMX is impractical for large systems.3.2.2 A practical preconditionerA good approximation for M +DTL−1D is required to make the ideal pre-conditioner, MiMX, suitable in practice. It has been shown in [27] thatM +DTL−1D is spectrally equivalent to M +X where X is the vector massmatrix on the magnetic space defined in (3.1). Using this approximationyields the practical preconditionerMMX =N 00 L, (3.8)where N = M + X is a shifted curl-curl operator. A scalable multigridsolver for N has been developed in [31] which involves one shifted Laplaciansolves on the vector space and one scalar Laplacian solve. However, the35construction of this multigrid solver is involved and hence we will considerdirect solves for this preconditioner and leave the multigrid implementationas a possible area of future work.3.3 Preconditioning the MHD equationsIn Section 2.3 and 2.4 we introduced three iteration schemes, namely Picarditeration (P), Magnetic Decoupling (MD) and Complete Decoupling (CD).Using the results from Sections 3.2 and 3.1 we will discuss the precondition-ing approaches that we apply for these non-linear iteration schemes.3.3.1 Picard iterationUsing the techniques of Sections 3.1 and 3.2 for the incompressible Navier-Stokes and Maxwell problems, respectively, we will look at possible scalablepreconditioners for the linearised MHD problemKMH =A+O BT CT 0B 0 0 0−C 0 M DT0 0 D 0. (3.9)36We propose the following preconditioner for KMHMMH =F BT CT 00 −S 0 0−C 0 N 00 0 0 L, (3.10)where S is now either the LSC or PCD approximation to the fluid flow Schurcomplement. The preconditioned matrix, M−1MHKMH, has eigenvalues λ = 1with algebraic multiplicity of at least nu + nb and the eigenvalue λ = −1with algebraic multiplicity of at least mb see [35, Theorem 8]. Due to thecoupling terms, C, the application of this preconditioner is computationallyexpensive. To overcome this, we propose to invert MMH by means of aninner preconditioner. The inner preconditioner is taken asMinnerMH =F BT 0 00 −S 0 00 0 N 00 0 0 L. (3.11)Here the preconditioned matrix, M−1innerMHMMH, has an eigenvalue λ = 1with algebraic multiplicity of at least mu + nu + 3mb − nb, see [35, Theo-rem 10].373.3.2 Magnetic DecouplingFrom Section 2.4.1 the matrix to be preconditioned for (MD) isKMD =F BT 0 0B 0 0 00 0 M DT0 0 D 0. (3.12)Recall that removing the coupling terms completely decouples the system.This therefore enables us to use the preconditioners for each of the subprob-lems separately and in parallel. Using the subproblem preconditioners (3.3)and (3.8) then we propose the following preconditioner for KMD:MMD =F BT 0 00 −S 0 00 0 N 00 0 0 L, (3.13)where S is the LSC or PCD approximation and N is the shifted curl-curlmatrix.383.3.3 Complete DecouplingFrom Section 2.4.2 the matrix to be preconditioned for (CD) isKCD =A BT 0 0B 0 0 00 0 M DT0 0 D 0. (3.14)Note that the matrix KCD is now symmetric. First we consider how to dealwith the upper (2, 2) block matrix which corresponds to the discrete StokesequationsKS =A BTB 0.As with the incompressible Navier-Stokes subproblem the idea for the Stokespreconditioner is again to approximate the Schur complementSS = BA−1BT .Recall that the matrix A is defined with the viscosity ν in Section 2.1. It wasshown in [56, 61] that the scaled pressure mass matrix, 1νW defined in (3.1),is spectrally equivalent to the Schur complement SS (which is also a conse-quence of the inf-sup stability from our mixed discretisation). Therefore apossible scalable Stokes preconditioner isA 00 1νW. (3.15)39Using (3.15) together with the Maxwell subproblem preconditioner (3.8)gives the preconditionerMCD =A 0 0 00 1νW 0 00 0 N 00 0 0 L. (3.16)As the matrix system KCD is symmetric, then the appropriate choice forthe Krylov subspace method is MINRES for each subproblem. The mainadvantage of using MINRES over GMRES is that we do not need to store anew vector at each iteration. Therefore, in terms of computational memoryMINRES is more efficient.3.3.4 SummaryIn summary, we outlined the three preconditioning approaches for the lin-earised systems arising in the non-linear iteration schemes proposed in Chap-ter 2. Table 3.1 references the coefficient matrices together the associatedpreconditioner.Iteration Coefficient (Outer)scheme matrix preconditioner(P) KMH in (3.9) MMH in (3.10)(MD) KMD in (3.12) MMD in (3.13)(CD) KCD in (3.14) MCD in (3.16)Table 3.1: Summary of coefficient matrices and corresponding precondition-ers for each iteration scheme40Note that for the Picard iteration (P), we employ the inner preconditionerMinnerMH in (3.11) to solve systems corresponding to (3.10).41Chapter 4Numerical resultsIn this chapter we present a series of convergence and preconditioning ex-periments. The principal aim is to check the scalability performance of thepreconditioned iterative methods for the MHD problem and the two decou-pling schemes proposed in Chapter 3. We have chosen to work with thePython programming language, due to its flexibility and the availability of alarge collection of external libraries that this environment offers. Python isan extremely portable platform, which allows for a nearly seamless interfacebetween multiple languages.4.1 SoftwareDue to the complex nature of the MHD problem, a number of differentlibraries have been used both to discretise and then to solve the resultingsystems.The finite element software that was used to discretise (1.11), (1.12)is FEniCS [40]. The core libraries used within FEniCS are the problem-solving interface DOLFIN [42, 43], the compiler for finite element variationalforms FFC [34, 41, 48], the finite element tabulator FIAT [32, 33] for creatingfinite element function spaces, the just-in-time compiler Instant, the code42generator UFC [2, 4] and the form language UFL [1, 3].Along with FEniCS we have used a number of stand-alone linear algebralibraries. The main one that has been used is PETSc4PY, which is the Pythoninterface for PETSc [9, 10]; it has been used mainly for the iterative solversas well as the blockwise preconditioning setup. The following packages wereused to solve the preconditioned system: HYPRE [19] as a multigrid solverand the sparse direct solvers UMFPACK [15, 16, 57, 58], PASTIX [29], SuperLU[37, 38] and MUMPS [5–7].4.2 Problem setupWhen numerically solving a set of linear or non-linear equations the problemneeds to be initialised: the mesh sequence, non-linear iteration stoppingcriteria and initial guess tolerance need to be defined. In this section, wewill briefly cover these aspects. For simplicity, all our numerical experimentsare based on lowest-order elements (i.e., k = 2).Mesh sequencesIn our tests, we consider sequences of uniformly refined simplicial grids(i.e., triangles in 2D and tetrahedra in 3D). The mesh level ` defines thenumber of edges between the nodes along a single boundary edge to be 2`.For example, consider the domain to be the unit square then the third level(` = 3) generates the uniform 8× 8 triangular grid given in Figure 4.1. Thefirst column of all convergence and iteration tables below will show the gridlevel.43Figure 4.1: Level 3 grid on unit square domainStopping criteriaThe incompressible Navier-Stokes equations as well as the full MHDproblem are a set of non-linear equations. Section 2.3 outlines the processin which we linearise the problem and then solve for the updates. Thestopping criteria that we enforce for the Navier-Stokes and MHD problemsare:• NS: ‖δu‖2 + ‖δp‖2 < tol,• MHD: ‖δu‖2 + ‖δp‖2 + ‖δb‖2 + ‖δr‖2 < tol,where ‖ · ‖2 denotes the 2-norm of a vector. The tolerance is stated in eachexample table.Note that the non-linear stopping tolerance that has been used variesfrom example to example. For the convergence results we set the tolerance44to be stricter to make sure the convergence rates are accurate. For the pre-conditioning tests the tolerance is loosened slightly for ease of computation.Initial guess toleranceTo initialise the non-linear iterations, we follow Section 2.5. As discussedthere, this involves iteratively solving Stokes equations (2.15) and a modifiedMaxwell problem of the form (2.16) to a sufficient accuracy. To that end,and throughout the chapter we choose the Krylov tolerance as 1e-10 for eachproblem. Inhomogeneous boundary conditions are enforced in a standardfashion within the corresponding finite element spaces, by interpolation atthe boundary degrees of freedom.The notation “Dofs” in the tables to follow stand for the number ofdegrees of freedom. We usually indicate which unknowns are meant butif not then “Dofs” refers to the total number of degrees of freedom of thesystem.4.3 Numerical results: Navier-Stokes equationsBefore considering the discretised MHD problem, we require that the Navier-Stokes subproblem perform as expected in terms of the error estimates andpreconditioner scalability. To check the validity of the FEniCS code, weintroduce two problems with smooth solutions, the first in 2D and the secondin 3D. As mention in Section 4.2, we only provide results lowest-order(k = 2) Taylor-Hood elements; the discretisation in this case is stable andis sufficiently accurate for our purposes.454.3.1 Convergence results for a 2D smooth solutionConsider the domain Ω = (0, 1)2. Taking the viscosity to be ν = 1, we choosethe source terms f and inhomogeneous Dirichlet boundary conditions on ∂Ωfrom the analytical solutionu(x, y) =sin(x) exp(x+ y) + cos(y) exp(x+ y)− sin(y) exp(x+ y),p(x, y) = x3 sin(y) + exp(x+ y).` Dofs uh/ph ‖u− uh‖L2(Ω) order ‖u− uh‖H1(Ω) order3 578/81 2.2100e-04 - 1.5645e-02 -4 2,178/289 2.7598e-05 3.00 3.9133e-03 2.005 8,450/1,089 3.4465e-06 3.00 9.7822e-04 2.006 33,282/4,225 4.3072e-07 3.00 2.4455e-04 2.007 132,098/16,641 5.3836e-08 3.00 6.1137e-05 2.008 526,338/66,049 6.7291e-09 3.00 1.5284e-05 2.00Table 4.1: Convergence of the velocity field for 2D Navier-Stokes - tol = 1e-10` Dofs uh/ph ‖p− ph‖L2(Ω) order3 578/81 7.2758e-03 -4 2,178/289 1.8190e-03 2.005 8,450/1,089 4.5025e-04 2.016 33,282/4,225 1.1341e-04 1.997 132,098/16,641 2.8639e-05 1.998 526,338/66,049 7.0146e-06 2.03Table 4.2: Convergence of the pressure variable for 2D Navier-Stokes -tol = 1e-10For the lowest-order Taylor-Hood elements the expected convergence46rates for the velocity field in the L2- and H1-norm errors are third andsecond order, respectively, and for the pressure the L2-norm error is secondorder. This is in agreement with observed orders from Tables 4.1 and 4.2.4.3.2 Convergence results for a 3D smooth solutionThe three-dimensional set up is very similar as with the 2D case. Thedomain is the unit cube Ω = (0, 1)3. As before, the viscosity is ν = 1, thenthe source term f and inhomogeneous Dirichlet boundary conditions on ∂Ωare calculated from the analytical solutionu(x, y, z) =− exp(x+ y + z) sin(y) + exp(x+ y + z) sin(z)exp(x+ y + z) sin(x)− exp(x+ y + z) sin(z)− exp(x+ y + z) sin(x) + exp(x+ y + z) sin(y),p(x, y, z) = exp(x+ y + z) + sin(y).Running the 3D code produces the results in Tables 4.3 and 4.4, which showsthe computed errors and convergence rates.` Dofs uh/ph ‖u− uh‖L2(Ω) order ‖u− uh‖H1(Ω) order1 375/27 2.6211e-02 - 4.5804e-01 -2 2,187/125 3.2997e-03 2.99 1.1547e-01 1.993 14,739/729 4.1267e-04 3.00 2.8944e-02 2.004 107,811/4,913 5.1565e-05 3.00 7.2416e-03 2.005 823,875/35,937 6.4443e-06 3.00 1.8108e-03 2.00Table 4.3: Convergence of the velocity field for 3D Navier-Stokes -tol = 1e105We see that the convergence rates are the same as for the 2D test solution47` Dofs uh/ph ‖p− ph‖L2(Ω) order1 375/27 1.7725e+00 -2 2,187/125 2.8602e-01 2.633 14,739/729 4.0587e-02 2.824 107,811/4,913 6.4794e-03 2.655 823,875/35,937 1.2724e-03 2.35Table 4.4: Convergence of the pressure variable for 3D Navier-Stokes -tol = 1e-10in the previous section for the velocity field. However, we observe that forthe pressure we get a slightly higher than expected order, namely about 2.5.4.3.3 Preconditioning with LSC and PCDIn this section, we check the performance of both the PCD and LSC pre-conditioners for the Navier-Stokes subproblem, see Sections 3.1.1 and 3.1.2.The problem setup is the 2D one from Section 4.3.1. We use GMRES asthe Krylov subspace solver (with a convergence tolerance of 1e-5) and theapplication of the preconditioners will be done directly. Note that, we donot restart GMRES. The results tables show GMRES iterations averagedover all non-linear iterations.We first consider the performance of the LSC preconditioner. In Table 4.5we show the average number of GMRES iterations for various values of ν andgrid levels. We observe that as the grid level increases the average numberof GMRES iterations increases. For the smaller mesh levels (` ≤ 4) theiteration increase is fairly small. However, for the higher levels the averageGMRES iterations seems to roughly double as the level doubles. This isbehaviour is similar to applying an incomplete factorisation preconditioner48` Dofs Average iterationsuh/ph ν = 10 ν = 1 ν = 0.1 ν = 0.011 50/9 5 5 6 82 162/25 10 10 10 213 578/81 16 16 15 304 2,178/289 25 24 24 315 8,450/1,089 48 47 46 396 33,282/4,225 103 86 77 637 132,098/16,641 190 190 141 131Table 4.5: Iteration table for LSC preconditioner for 2D example for variousvalues of ν and tol = 1e-5to the problem. Hence, LSC does not seem to yield a scalable preconditioner.We next look at how PCD performs as a preconditioner for the Navier-Stokes subproblem. Applying this preconditioner to the incompressibleNavier-Stokes problem in isolation gives the results in Table 4.6. We ob-serve that as the mesh level increases, the average number of GMRES itera-tion numbers stays roughly constant with respect to the mesh level `. Thisscalability is what we expect and require for the Navier-Stokes subproblempreconditioner. Also, note for fixed values of ν the number of average num-ber of GMRES iterations remains about constant. We see that for ` = 5, 6and ν = 0.01 the iterations seem to be slightly higher than expected. Thisbehaviour is likely due to the mesh size being too big.4.4 Numerical results: Maxwell’s equationsIn this section we consider 2D and 3D test problems to validate the FEniCScode for the Maxwell subproblem in mixed form.49` Dofs Average iterationsuh/ph ν = 10 ν = 1 ν = 0.1 ν = 0.015 8,450/1,089 17 17 21 586 33,282/4,225 17 17 22 307 132,098/16,641 18 17 22 218 526,338/66,049 18 18 22 209 2,101,250/263,169 18 19 22 21Table 4.6: Iteration table for a PCD preconditioned 2D example for variousvalues of ν and tol = 1e-54.4.1 Convergence results for a 2D smooth solutionAgain, consider the domain Ω = (0, 1)2 with κ = νm = 1. We choose thesource term g and inhomogeneous Dirichlet boundary conditions on ∂Ω fromthe analytical solutionb(x, y) =exp(x+ y) cos(x)exp(x+ y) sin(x)− exp(x+ y) cos(x),r(x, y) = sin(2pix) sin(2piy).` Dofs bh/rh ‖b− bh‖L2(Ω) order ‖b− bh‖H(curl,Ω) order1 48/25 9.3833e-02 - 1.0696e-01 -2 176/81 2.3350e-02 2.01 2.7088e-02 1.983 672/289 5.8372e-03 2.00 6.4220e-03 2.084 2,624/1,089 1.4597e-03 2.00 1.4586e-03 2.145 10,368/4,225 3.6497e-04 2.00 3.5066e-04 2.066 41,216/16,641 9.1246e-05 2.00 8.6688e-05 2.027 164,352/66,049 2.2812e-05 2.00 2.1609e-05 2.00Table 4.7: Convergence for 2D Maxwell smooth solution - magnetic fieldFor Maxwell’s equations in mixed form we use second order (k = 2)50` Dofs bh/rh ‖r − rh‖L2(Ω) order ‖r − rh‖H1(Ω) order1 48/25 2.7761e-01 - 2.8932e+00 -2 176/81 4.3540e-02 2.67 9.3299e-01 1.633 672/289 4.8633e-03 3.16 2.5904e-01 1.854 2,624/1,089 5.6724e-04 3.10 6.6810e-02 1.965 10,368/4,225 6.9363e-05 3.03 1.6841e-02 1.996 41,216/16,641 8.6203e-06 3.01 4.2193e-03 2.007 164,352/66,049 1.0760e-06 3.00 1.0554e-03 2.00Table 4.8: Convergence for 2D Maxwell smooth solution - multiplier variableNe´de´lec elements of the first kind [47] and quadratic nodal elements for themagnetic field and multiplier, respectively. We expect second order conver-gence in the errors ‖b − bh‖L2(Ω), ‖b − bh‖H(curl,Ω) and ‖r − rh‖H1(Ω), andthird order in ‖r − rh‖L2(Ω), see [44]. Tables 4.7 and 4.8 confirm that theserates are indeed obtained for both the magnetic field and the multiplier.4.4.2 Convergence results for a 3D smooth solutionAs we did for the 3D Navier-Stokes problem the domain is the unit cube,i.e. Ω = (0, 1)3, κ = νm = 1 and using the analytical solutionb(x, y, z) =− exp(x+ y + z) sin(y) + exp(x+ y + z) sin(z)exp(x+ y + z) sin(x)− exp(x+ y + z) sin(z)− exp(x+ y + z) sin(x) + exp(x+ y + z) sin(y),r(x, y, z) = sin(2pix) sin(2piy) sin(2piz),the source term g and inhomogeneous Dirichlet boundary conditions on ∂Ωare defined.Tables 4.9 and 4.10 show that the orders of convergence do asymptoti-51` Dofs bh/rh ‖b− bh‖L2(Ω) order ‖b− bh‖H(curl,Ω) order1 436/125 5.3312e-02 - 3.5258e-01 -2 2,936/729 1.4192e-02 1.91 8.9944e-02 1.973 21,424/4,913 3.5801e-03 1.99 2.2690e-02 1.994 163,424/35,937 8.9697e-04 2.00 5.7585e-03 1.98Table 4.9: Convergence for 3D Maxwell smooth solution - magnetic field` Dofs bh/rh ‖r − rh‖L2(Ω) order ‖r − rh‖H1(Ω) order1 436/125 2.2689e-01 - 2.8923e+00 -2 2,936/729 5.9452e-02 1.93 1.1782e+00 1.303 21,424/4,913 6.9068e-03 3.11 3.3901e-01 1.804 163,424/35,937 7.5913e-04 3.19 9.0082e-02 1.91Table 4.10: Convergence for 3D Maxwell smooth solution - multiplier vari-ablecally approach the expected orders as in the 2D case from Section 4.4.1.4.4.3 PreconditioningSection 3.2 introduces the preconditioning strategy for Maxwell’s equationsthat will be employed in this thesis. Two- and three-dimensional results arepresented in this section.2D exampleTable 4.11 shows the performance of using the preconditioner MMX definedin (3.8) with a MINRES tolerance of 1e-6 and a direct application of thepreconditioner for the 2D problem in Section 4.4.1.From Table 4.11 we observe that the number of MINRES iterations staysabout constant as the mesh level, `, increases. Note that for a fixed magneticviscosity (νm), the number of iterations remains roughly the same, and the52` Dofs bh/rh Number of iterationsνm = 10 νm = 100 νm = 1000 νm = 100005 10,368/4,225 5 4 6 66 41,216/16,641 5 6 6 67 164,352/66,049 5 6 6 68 656,384/263,169 5 6 6 89 2,623,488/1,050,625 4 6 6 810 10,489,856/4,198,401 4 6 8 10Table 4.11: Iteration count for Maxwell preconditioner for 2D example -direct application of preconditionerperformance is fairly insensitive with respect to the values of νm.3D exampleRunning the 3D code for the problem in Section 4.4.2 produces the results inTable 4.12. Again, the number of MINRES iterations stay roughly constantas the mesh level, `, increases. We also note that the number of iterationsseems again to be insensitive to the value of the magnetic viscosity, νm.` Dofs bh/rh Number of iterationsνm = 10 νm = 100 νm = 1000 νm = 100002 2,936/729 4 4 6 53 21,424/49,13 4 4 6 54 163,424/35,937 4 6 6 55 1,276,096/274,625 4 6 6 5Table 4.12: Iteration table for Maxwell preconditioner for 3D example -direct application of preconditioner534.5 Numerical results: MHD problemSections 4.3 and 4.4 show results for the Navier-Stokes and Maxwell’s equa-tions in isolation. The next step is to incorporate these two subproblemsinto the full MHD system.4.5.1 Convergence results for a 2D smooth solutionTo validate the code, we consider the following 2D test problem: let thedomain be the unit square Ω = (0, 1)2. Take ν = κ = 1, νm = 10; then thesource terms f , g and inhomogeneous Dirichlet boundary conditions on ∂Ωare defined from the analytical solution:u(x, y) =xy exp(x+ y) + x exp(x+ y)−xy exp(x+ y)− y exp(x+ y),p(x, y) = exp(y) sin(x),b(x, y) =exp(x+ y) cos(x)exp(x+ y) sin(x)− exp(x+ y) cos(x),r(x, y) = x sin(2pix) sin(2piy).The convergence results computed in Tables 4.13 to 4.16 agree with theoptimal rates for the errors ‖u − uh‖L2(Ω), ‖u − uh‖H1(Ω), ‖b − bh‖L2(Ω),‖b−bh‖H(curl,Ω), ‖r−rh‖L2(Ω) and ‖r−rh‖H1(Ω). The pressure convergencerate exhibits a little higher than expected rate of convergence for ` < 5 buteventually settles down to second order for the higher levels.54` Dofs uh/ph ‖u− uh‖L2(Ω) order ‖u− uh‖H1(Ω) order3 578/81 1.3357e-03 - 8.3244e-02 -4 2,178/289 1.6752e-04 3.00 2.0840e-02 2.005 8,450/1,089 2.1029e-05 2.99 5.2110e-03 2.006 33,282/4,225 2.6675e-06 2.98 1.3028e-03 2.007 132,098/16,641 3.5200e-07 2.92 3.2570e-04 2.008 526,338/66,049 5.2267e-08 2.75 8.1425e-05 2.00Table 4.13: Convergence for 2D MHD smooth solution (velocity field) -tol = 1e-8` Dofs uh/ph ‖p− ph‖L2(Ω) order3 578/81 4.9329e-03 -4 2,178/289 7.4197e-04 2.735 8,450/1,089 1.5604e-04 2.256 33,282/4,225 1.5604e-04 2.257 132,098/16,641 3.7738e-05 2.058 526,338/66,049 2.3415e-06 2.00Table 4.14: Convergence for 2D MHD smooth solution (pressure variable) -tol = 1e-8` Dofs bh/rh ‖b− bh‖L2(Ω) order ‖b− bh‖H(curl,Ω) order3 672/289 7.1735e-03 - 1.1521e-02 -4 2,624/1,089 1.7944e-03 2.00 2.8833e-03 2.005 10,368/4,225 4.4868e-04 2.00 7.2102e-04 2.006 41,216/16,641 1.1217e-04 2.00 1.8027e-04 2.007 164,352/66,049 2.8044e-05 2.00 4.5067e-05 2.008 656,384/263,169 7.0110e-06 2.00 1.1267e-05 2.00Table 4.15: Convergence for 2D MHD smooth solution (magnetic field) -tol = 1e-855` Dofs bh/rh ‖r − rh‖L2(Ω) order ‖r − rh‖H1(Ω) order3 672/289 3.1956e-03 - 1.6711e-01 -4 2,624/1,089 3.6326e-04 3.14 4.2565e-02 1.975 10,368/4,225 4.3989e-05 3.05 1.0692e-02 1.996 41,216/16,641 5.4517e-06 3.01 2.6763e-03 2.007 164,352/66,049 6.7998e-07 3.00 6.6928e-04 2.008 656,384/263,169 8.4951e-08 3.00 1.6733e-04 2.00Table 4.16: Convergence for 2D MHD smooth solution (multiplier variable)- tol = 1e-84.5.2 Convergence results for a 3D smooth solutionNewt, we consider a 3D test problem discretised on the unit cube Ω = (0, 1)3with Dirichlet boundary conditions. Let ν = κ = 1, νm = 10 and let theanalytical solution be given by:u(x, y, z) =−xy exp(x+ y + z) + xz exp(x+ y + z))xy exp(x+ y + z)− yz exp(x+ y + z)−xz exp(x+ y + z) + yz exp(x+ y + z),p(x, y, z) = exp(x+ y + z) sin(y),b(x, y, z) =− exp(x+ y + z) sin(y) + exp(x+ y + z) sin(z)xy exp(x+ y + z)− yz exp(x+ y + z)− exp(x+ y + z) sin(x) + exp(x+ y + z) sin(y),r(x, y, z) = sin(2pix) sin(2piy) sin(2piz).Then the source terms f and g and inhomogeneous boundary conditions aredefined from the analytical solution.From Tables 4.17 to 4.20 it can be seen that we obtain roughly the same56` Dofs uh/ph ‖u− uh‖L2(Ω) order ‖u− uh‖H1(Ω) order1 375/27 4.2196e-03 - 7.8282e-02 -2 2,187/125 5.2733e-04 3.00 1.9495e-02 2.013 14,739/729 6.5749e-05 3.00 4.8664e-03 2.004 107,811/4,913 8.2092e-06 3.00 1.2161e-03 2.00Table 4.17: Convergence for 3D MHD smooth solution (velocity field) -tol = 1e-8` Dofs uh/ph ‖p− ph‖L2(Ω) order1 375/27 5.0614e-02 -2 2,187/125 9.0306e-03 2.493 14,739/729 1.9035e-03 2.254 107,811/4,913 4.5311e-04 2.07Table 4.18: Convergence for 3D MHD smooth solution (pressure variable) -tol = 1e-8` Dofs bh/rh ‖b− bh‖L2(Ω) order ‖b− bh‖H(curl,Ω) order1 436/125 3.9931e-03 - 2.4817e-02 -2 2,936/729 1.1211e-03 1.83 7.4101e-03 1.743 21,424/4,913 2.9642e-04 1.92 1.9439e-03 1.934 163,424/35,937 7.5494e-05 1.97 4.9180e-04 1.98Table 4.19: Convergence for 3D MHD smooth solution (magnetic field) -tol = 1e-8` Dofs bh/rh ‖r − rh‖L2(Ω) order ‖r − rh‖H1(Ω) order1 436/125 7.4750e-04 - 1.0113e-02 -2 2,936/729 9.2457e-05 3.02 2.9368e-03 1.783 21,424/4,913 1.1182e-05 3.05 7.7197e-04 1.934 163,424/35,937 1.3829e-06 3.02 1.9597e-04 1.98Table 4.20: Convergence for 3D MDH smooth solution (multiplier variable)- tol = 1e-857convergence orders as with the 2D case. The computed orders are not quiteas clean as with the 2D results and in some cases the calculated order isonly slightly lower than the expected rate. However, for mesh levels 3 and4 the computed‘’ orders seem to converge to the expected rates.4.5.3 Parameter testsThere are two different decoupling schemes considered in thesis ((MD) and(CD)). The convergence of the non-linear iterations (within some tolerancetol) is likely to be limited by the parameter setup of the problem, i.e., by thevalues of the fluid viscosity (ν), the magnetic viscosity (νm) and the couplingnumber (κ). In this section we will look at how these two schemes performwith respect to the Picard iteration (P) when fixing two parameters andvarying the third. We will only look at varying κ and ν since νm appearswithin all three schemes. Also, if the non-linear iterations do not convergethen this is denoted by “-” in the tables.Viscosity testAs a first test we consider varying the fluid viscosity, ν, for tol = 1e-5, κ = 1and νm = 10. The results are shown in Table 4.21.As the fluid viscosity (ν) decreases, the fluid flow equations (1.11a) and(1.11b) become more convection-dominated. Thus we see that the (CD)scheme breaks down for small ν, as in this decoupling scheme the convectionterm is taken explicitly. On the other hand, the Picard and (MD) schemesperform similarly. For ν = 0.0001 all schemes break down and do notconverge. This may be due to the stability of the convection discretisation58ν = 1 ν = 0.1 ν = 0.01 ν = 0.001` Dofs (P) (MD) (CD) (P) (MD) (CD) (P) (MD) (CD) (P) (MD) (CD)3 1,620 5 7 15 8 10 - 12 14 - 40 40 -4 6,180 5 7 15 8 10 - 12 14 - 40 40 -5 24,132 5 7 15 8 10 - 12 14 - 19 19 -6 95,364 5 7 15 8 10 - 12 14 - 17 19 -7 379,140 5 7 15 8 10 - 12 14 - 17 19 -8 1,511,940 5 7 15 8 10 - 12 14 - 17 19 -Table 4.21: Number of non-linear iterations for various values of ν withtol = 1e-5, κ = 1 and νm = 10.or it may be that the solution is unstable for these Reynolds numbers.Coupling number testκ = 1 κ = 10 κ = 100 κ = 1000` Dofs (P) (MD) (CD) (P) (MD) (CD) (P) (MD) (CD) (P) (MD) (CD)3 1,620 5 5 13 6 7 17 7 13 26 7 93 -4 6,180 5 5 15 6 7 18 7 13 26 7 93 -5 24,132 5 5 15 6 7 18 7 13 26 7 93 -6 95,364 5 5 15 6 7 18 7 13 26 7 93 -7 379,140 5 5 15 6 7 18 7 13 26 7 93 -8 1,511,940 5 5 15 6 7 18 7 13 26 7 93 -Table 4.22: Number of non-linear iterations for various values of κ withtol = 1e-5, ν = 1 and νm = 10.The next parameter test examines the effects of the coupling terms inthe three non-linear iteration schemes. We expect the Picard scheme tooutperform the other schemes for large values of κ. The results in Table 4.22show that this is indeed the case. For the 2D example in Section 4.5.1 withν = 1, νm = 10, the (CD) scheme completely breaks down for κ > 100,whereas the number of non-linear iterations for (MD) stays roughly constantfor each mesh level but is larger than for the Picard iteration. When κ =591000, (MD) also converges more slowly. This is the point at which the Picarditeration (P) becomes the most viable option.4.6 Preconditioned MHD problemSection 3.3 introduced the three preconditioning techniques that we applyto each of the decoupling schemes and the Picard iteration. In this sectionwe investigate in detail the performance of these three strategies.In the subsequent tables we use the following notation:• ItsNL: the number of non-linear iterations (used in (P), (MD) and(CD));• ItsO: the average number of outer GMRES iterations (used in (P));• ItsI: the average number of inner iterations (used in (P));• ItsNS: the average number of iterations to solve the Navier-Stokesproblem (used in (MD));• ItsS: the average number of iterations to solve the Stokes problem(used in (CD));• ItsM: the average number of iterations to solve the Maxwell problem(used in (MD) and (CD));• Dofs: the total degrees of freedom of the system.• Av solve time: the solve time (with out matrix assembly) averagedover all non-linear iterations60• Total time: the total solve and assembly time until the non-lineariterations convergeThe problem setups are the same as Sections 4.5.1 and 4.5.2. The numericaltests presented in this section are two-dimensional unless otherwise stated.As mentioned in Chapter 3, the application of the preconditioners is be donethrough direct solves.4.6.1 Picard iterationIn this section, we present preconditioning results for the Picard iteration(P) given in Section 3.3.1. Recall, this is the inner-outer preconditioningtechnique. The only difference between the outer preconditioner (3.10)and inner preconditioner (3.11) is the inclusion of the coupling terms Cand −CT . Therefore, it may be possible to relax the inner solve toleranceand still produce scalable results with a smaller number of inner solves. Intougher cases (e.g., large κ) we expect to have a trade off, that is, an innertolerance that is too loose may degrade convergence. We start off withparameter tests similar to those in Section 4.5.3.Preconditioned parameter testsThe first preconditioning parameter test we consider is again to see how theviscosity, ν, effects the preconditioners performance for κ = 1 and νm = 10.For the preconditioner for the Picard iteration we apply the inner-outerpreconditioning approach from Section 3.3.1. For the outer Krylov solverwe use flexible GMRES (FGMRES), see [52], with a convergence tolerance61of 1e-6 and for the inner Krylov solver we use GMRES with a tolerance of1e-6. The results are presented in Table 4.23.ν = 0.1 ν = 1 ν = 10` Dofs ItsNL ItsO ItsI ItsNL ItsO ItsI ItsNL ItsO ItsI4 6,180 - - - 6 18.8 14.7 5 11.2 12.65 24,132 - - - 6 27.3 15.6 6 9.8 12.76 95,364 13 91.5 63.4 6 20.2 15.3 6 9.7 11.37 379,140 10 74.6 31.4 7 16.3 14.3 7 14.6 15.98 1,511,940 12 82.3 37.8 8 24.3 14.9 7 15.4 15.7Table 4.23: Number of non-linear and average number of preconditioningiterations for various values of ν with tol = 1e-5, κ = 1 and νm = 10.We observe that for smaller values of ν (i.e., in the convection-dominatedregime) the number of inner and outer iterations increases. The table showsthat as the mesh level increases the number of inner and outer iterations stayroughly constant. We also note that for mesh levels 4 and 5 the iterationsdo not seem to converge for ν = 0.1. This may be because we are notconsidering grid levels of sufficient size. The PCD approximation to theSchur complement becomes more accurate for when h decreases.The second test we consider is to again vary the coupling number. Asbefore, we use the outer convergence tolerance of 1e-6 and inner tolerance1e-6 with the same Krylov solvers.Table 4.24 shows that as the mesh level increases then the number ofouter and inner GMRES iterations stays roughly constant. We note thatfor larger values of κ the number of outer iterations increases. We also seethis trend with the inner iterations for κ = 10 and 100.62κ = 0.1 κ = 1 κ = 10 κ = 100` Dofs ItsNL ItsO ItsI ItsNL ItsO ItsI ItsNL ItsO ItsI ItsNL ItsO ItsI4 6,180 5 22.2 15.4 6 18.8 14.7 7 30.7 14.1 9 61.4 24.45 24,132 5 23.8 11.4 6 27.3 15.2 8 43.8 24.0 10 80.3 37.96 95,364 5 28.0 17.4 6 20.2 15.3 7 41.1 15.4 9 74.6 31.17 379,140 5 18.6 14.6 7 16.3 14.2 7 37.4 16.4 14 73.9 34.78 1,511,940 5 20.4 15.2 8 24.3 14.9 7 39.6 18.4 11 75.4 33.3Table 4.24: Number of non-linear and average number of preconditioningiterations for various values of κ with tol = 1e-5, ν = 1 and νm = 10.4.6.2 Magnetic DecouplingThis section will look at preconditioning results for the (MD) scheme out-lined in Section 3.3.2. As with the Picard iterations results we will considerparameter tests to assess the robustness of this decoupling scheme. Alltests considered here will have a GMRES stopping tolerance of 1e-6 for theNavier-Stokes subproblem and a MINRES stopping tolerance of 1e-6 for theMaxwell subproblem.Preconditioned parameter testsWe will again start of by seeing how the fluid’s viscosity affects the precon-ditioner. The results are shown in Table 4.25.ν = 0.01 ν = 0.1 ν = 1 ν = 10` Dofs ItsNL ItsNS ItsM ItsNL ItsNS ItsM ItsNL ItsNS ItsM ItsNL ItsNS ItsM5 24,132 13 360.3 3.3 8 24.1 3.3 5 22.0 3.4 4 22.0 3.56 95,364 13 105.2 3.3 8 25.1 3.4 5 21.2 3.4 4 20.8 3.57 379,140 13 71.5 3.4 8 25.3 3.3 5 21.2 3.4 4 21.0 3.58 1,511,940 13 65.3 3.4 8 25.9 3.4 5 21.4 3.4 4 21.3 3.5Table 4.25: Number of non-linear iterations and average number of iterationsto solve the Navier-Stokes and Maxwell’s subproblem for the MD schemewith tol = 1e-4, κ = 1, ν = 1 and νm = 10.63From Table 4.25 we see that for ν = 0.01 the Navier-Stokes solver con-verges slowly to the solution for the lower mesh levels (` ≤ 6). A possiblereason for this is that the PCD preconditioner for the Navier-Stokes equa-tions does not always converge for small viscosities and small meshes (largemesh size h). However, we see that for ν = 0.1, ν = 1 and ν = 10 both theNavier-Stokes and Maxwell’s subproblems exhibit iterations independent ofthe mesh level.κ = 0.1 κ = 1 κ = 10 κ = 100` Dofs ItsNL ItsNS ItsM ItsNL ItsNS ItsM ItsNL ItsNS ItsM ItsNL ItsNS ItsM5 24,132 4 22.1 4.5 5 22.0 3.4 10 21.4 2.3 - 21.7 2.26 95,364 4 21.5 4.5 5 21.2 3.4 10 21.1 2.3 - 21.5 2.37 379,140 4 21.5 4.8 5 21.2 3.4 10 21.1 2.4 - 21.6 2.28 1,511,940 4 21.5 4.8 5 21.4 3.4 10 21.1 2.4 - 21.7 2.2Table 4.26: Number of non-linear iterations and average number of iterationsto solve the Navier-Stokes and Maxwell’s subproblem for the MD schemewith tol = 1e-4, ν = 1 and νm = 10.Table 4.26 shows the iteration count of the Magnetic Decoupling schemewhen the coupling number, κ, is varied. We see that as κ increases thenumber of non-linear iterations increases. In fact, for κ = 100 the schemedoes not converge. If we recall Table 4.22 we note that for direct solvesκ = 100 does converge to the solution. Therefore, it seems that the Krylovtolerance of 1e-6 is not tight enough for (MD) to converge.Table 4.27 shows large scaled results for the (MD) scheme with up to24 million degrees of freedom. For all mesh levels the Navier-Stokes andMaxwell solves take roughly 21 and 3 iterations each. We also note that for` = 10, the total solve time was a little under 5 and a half hours.64` Dofs Av solve time Total time ItsNL ItsNS ItsM5 24,132 0.5 7.3 5 22.0 3.46 95,364 2.5 30.4 5 21.2 3.47 379,140 13.1 134.7 5 21.2 3.48 1,511,940 69.8 627.3 5 21.4 3.49 6,038,532 407.7 3159.7 5 21.6 3.210 24,135,684 3022.1 19668.3 5 21.6 3.4Table 4.27: Number of non-linear iterations and average number of iterationsto solve the Navier-Stokes and Maxwell’s subproblem for the MD schemewith tol = 1e-4, κ = 1, ν = 1 and νm = 10.Preconditioned 3D resultsTests were also run in three-dimensions and the results are represented be-low. Table 4.28 lists the number of iterations to solve the Navier-Stokesand Maxwell subproblems. The results seem to be independent of the meshlevel.l Dofs Av solve time Total time ItsNL ItsNS ItsM1 963 0.04 2.7 5 23.2 3.42 5977 0.20 17.2 5 34.2 3.03 41805 4.86 151.1 5 34.2 3.44 312,085 242.1 2222.8 5 32.4 3.45 2,410,533 30222.3 159032.8 5 30.8 3.2Table 4.28: Number of non-linear iterations and average number of iterationsto solve the Navier-Stokes and Maxwell’s subproblem for the MD schemewith tol = 1e-4, κ = 1, ν = 1 and νm = 10 in 3D.4.6.3 Complete DecouplingFinally, we consider preconditioning the Complete Decoupling scheme pre-sented in Section 3.3.3. As for (P) and (MD) we will consider varying boththe fluid’s viscosity as well as the coupling number. To solve the Stokes65block of (3.14) we will use MINRES, as the Stokes equations are symmetric.Results are shown with 1e-6 as the convergence tolerance for both the Stokesand Maxwell solvers for all tests.Preconditioned parameter testsTable 4.21 shows that the (CD) iteration scheme breaks down for ν = 0.1.Therefore, we consider tests for ν = 1 and ν = 10 and the results are givenin Table 4.29. From Table 4.29 we see that the number average iterations forν = 1 ν = 10` Dofs ItsNL ItsS ItsM ItsNL ItsS ItsM5 24,132 11 29.0 3.3 4 29.0 3.56 95,364 11 28.5 3.4 4 29.0 3.57 379,140 11 27.5 3.4 4 29.0 3.58 1,511,940 11 28.3 3.5 4 29.0 3.5Table 4.29: Number of non-linear iterations and average number of itera-tions to solve the Stokes and Maxwell’s subproblem for the CD scheme withtol = 1e-4 κ = 1, and νm = 10.the Stokes and Maxwell subproblem solves remains constant. This is to beexpected as we are in a regime (large enough ν) where the (CD) iterationsconverge. We also note that for larger values of ν the number of non-lineariterations increases.To finish the parameter tests in this thesis, we again consider varyingthe coupling number. The results in Table 4.30 shows that for κ = 0.1 to 10the number of iterations to solve the Stokes and Maxwell’s problems seemboth independent of the mesh size and the coupling number. The numberof non-linear iterations increase for larger κ, in fact, for κ = 100 the schemedoes not converge to the solution. We again consult with Table 4.22 where66κ = 0.1 κ = 1 κ = 10 κ = 100` Dofs ItsNL ItsS ItsM ItsNL ItsS ItsM ItsNL ItsS ItsM ItsNL ItsS ItsM5 24,132 9 29.0 4.4 11 29.0 3.3 18 29.0 2.2 - - -6 95,364 9 29.0 4.6 11 28.5 3.4 18 29.0 2.3 - - -7 379,140 8 28.5 4.5 11 27.6 3.4 18 29.0 2.3 - - -8 1,511,940 8 28.3 4.6 11 28.3 3.5 18 29.0 2.4 - - -Table 4.30: Number of non-linear iterations and average number of itera-tions to solve the Stokes and Maxwells subproblem for the CD scheme withtol = 1e-4, ν = 1 and νm = 10.a direct application of the coefficient matrices was considered and we seethat the corresponding iteration test there does converge for κ = 100. Wetherefore speculate, as for (MD), that a Krylov tolerance of 1e-6 is not tightenough when κ = 100.As with the (MD) scheme we test the (CD) iterations in a large scalemanner. From Table 4.31 we see that for l = 10 the total time is approx-imately 7 hours 50 minutes where as the for the same level using (MD) itwas 5 hours 30 minutes. This increase in time, even so the average solvetime is less, is due the fact that (CD) requires 6 more non-linear iterationsto converge to the solution.` Dofs Av solve time Total time ItsNL ItsS ItsM5 24,132 0.4 7.7 11 29.0 3.36 95,364 2.6 38.8 11 28.5 3.47 379,140 13.0 181.7 11 27.5 3.48 1,511,940 66.5 888.0 11 28.3 3.59 6,038,532 358.0 4565.4 11 28.3 3.410 24,135,684 2335.7 28337.4 11 27.4 3.5Table 4.31: Number of non-linear iterations and average number of itera-tions to solve the Stokes and Maxwell’s subproblem for the CD scheme withtol = 1e-4, κ = 1, ν = 1 and νm = 10.67Preconditioned 3D resultsThe final test run in this thesis is a three-dimensional example using the(CD) scheme. The results are shown in Table 4.32. This table, as with the2D results, show that the average iterations to solve Stokes and Maxwellsubproblems are independent with respect to the mesh level.` Dofs Av solve time Total time ItsNL ItsS ItsM1 963 0.03 1.8 6 30.0 3.52 5,977 0.20 9.9 6 45.7 3.23 41,805 4.32 89.2 6 43.0 2.84 312,085 214.3 1786.5 6 42.3 2.85 2,410,533 26954.4 165671.1 6 41.3 2.8Table 4.32: Number of non-linear iterations and average number of itera-tions to solve the Stokes and Maxwell’s subproblem for the CD scheme withtol = 1e-4, κ = 1, ν = 1 and νm = 10 in 3D.68Chapter 5Conclusions and future work5.1 ConclusionsIn this thesis, we developed an approach for the numerical solution of anincompressible magnetohydrodynamics (MHD) model, with the goal of solv-ing large scale problems. To that end, we generated a large scale code base,utilising both the finite element software package FEniCS [40] and linearalgebra software PETSc [9, 10].Our mixed finite element discretisation approach is based on usingTaylor-Hood elements for the fluid variables and on a mixed Ne´de´lec pairfor the magnetic unknowns. The three linearised iteration strategies for theMHD model range from standard Picard iterations to completely decoupledschemes. For these iteration schemes, we followed the preliminary precondi-tioning results from [35]. The proposed preconditioning ideas are motivatedby state of the art preconditioners for the mixed Maxwell and Navier-Stokessubproblems.We have introduced two preconditioners for the Navier-Stokes equations.When applying the Least-Squares Commutator preconditioner we observeda lack of scalability with respect to the mesh size and the viscosity parameterbut for the Pressure Convection-Diffusion (PCD) preconditioner we obtained69scalable results. This lead us to using PCD within the preconditioner forthe non-linear iteration schemes. The preconditioning results for Maxwell’sequations exhibited iteration counts independent of the mesh size and themagnetic viscosity.We have presented results for the three iteration schemes the linearisedMHD model. After confirming that the numerical results meet the expecteda priori convergences rates, parameter tests (using direct solvers for thematrix coefficients of the iteration schemes) were carried out to see the per-formance of the two decoupling schemes as well as the Picard iteration.From the parameter tests we saw that the robustness of the three schemesgreatly depends on the non-dimensional parameters in the problem. The Pi-card iterations behaved similarly to the Magnetic Decoupling (MD) schemefor relevant parameters based on physical applications. However, the Com-plete Decoupling (CD) scheme did not converge to the solution for smallkinematic viscosities (i.e., ν = 0.1).Finally, we showed preconditioning results for the MHD linearisationsresulting from the three non-linear schemes. One of the principal objectivesof this thesis was to perform large scale preconditioning tests. Here 2Dscalable experiments were run with 24 million degrees of freedom for the twodecoupled iterations. For the Picard iteration, results were only presentedfor up to 1.5 million degrees of freedom. We were able to go to muchgreater problem sizes for the decoupled iterations ((MD) and (CD)) as thepreconditioning approach is much less involved (no inner-outer iterations).We also considered the same parameter tests as above but using thepreconditioning approaches; the results examined the robustness of the de-70coupled preconditioned techniques. Again, the numerical results showediterations counts independent of the problem size. From the results theconvergence of all three iterations ((P), (MD) and (CD)) was hindered bysmall values of ν. In fact, only the (MD) scheme converged for a ν = 0.01.We stress the need to find more robust preconditioning techniques for smallvalues of the viscosity; see Section 5.2. Similar conclusions can be drawnfrom the scalable 3D numerical results presented here.From the preconditioned and non-preconditioned parameter tests in thisthesis, it seems that (MD) performs best for feasible physical parameters(small values of ν, κ ≈ 1 and νm ≈ 10). The Picard iteration is significantlymore costly due to the inner-outer iterations and hence, should only be usedin the ranges where (MD) or (CD) do not converge. If memory and thevalues of the parameters are moderate (e.g., ν ≥ 1) then the (CD) schemeis best. This is because the coefficient matrix for (CD) is symmetric andtherefore we use the cheaper MINRES Krylov solver.The experiments carried out in the thesis are quite idealised. We em-ployed uniform regular grids, and tested problems for smooth solutions withrelatively harmless parameter values. For many practical applications theproblems may involve complex 3D geometries with varying, discontinuousor even non-linear material coefficients. In addition, MHD flows may beconvection-dominated. Tackling these difficulties in a computationally effi-cient manner may require additional computational tools such as automaticmesh generation in complicated 3D domains, adaptive mesh refinement tech-niques, as well as refined finite element models to handle strong convectionand more complex material laws.71As for solution methods, in this thesis we have considered problems ofsufficiently small size which enabled us to use direct solves on the precon-ditioners. In practice, however, large real-world 3D problems require fullyscalable preconditioned iterative solution methods, e.g., specialised multi-grid techniquesThe scalable results presented in this thesis provide a starting point forfurther in-depth studies into these issues in the context of MHD applications.5.2 Future workWe finish this thesis with some considerations of possible areas for futurework.1. Scalable inner solvers:In this thesis, the application of the preconditioners have been carriedout using direct solves for the component blocks. However, to con-struct fully scalable preconditioners we require an iterative approachto solve these systems of equations. For the most part, this can bereadily achieved by using standard algebraic or geometric multigridcycles and/or iterative methods. We note though that the scalableiterative solution for the shifted curl-curl operator in (3.8) is more in-volved. One option is to employ the multigrid solver developed in [31].We expect the implementation of scalable inner solves to significantlyimprove the efficiency of the overall approach, especially for large scale3D problems and we should obtain a fully scalable iterative solver.722. Release code on a public repository:Once the code has been cleaned up and properly commented, the ideais to release it for public use on either GitHub https://github.com/or Bitbucket https://bitbucket.org/.3. Parallelisation of the code:The code to discretise and solve the MHD system (1.11), (1.12) iscurrently executed sequentially. For 3D problems of large dimension,an efficient parallel implementation is expected to greatly decreasesetup and solve time, in particular for the decoupling schemes.4. Robustness with respect to kinematic viscosity:The results presented in Chapter 4 indicate the efficient performancewith respect to the mesh size and with respect to the non-dimensionalparameters (ν, νm and κ). However, for small fluid viscosities, theperformance of all the proposed preconditioners for our discretisationstarts to degrade. The development of preconditioners and convectionstabilised discretisations that work well together and are more robustwith respect to ν is another area of possible future work. The solutionof convection-dominated flow is very challenging and is a subject ofactive research in many other areas.5. High frequencies in Maxwell’s equations:The Maxwell formulation considered in this MHD model only con-siders low frequencies. The consideration of high frequencies intro-duces hyperbolicity into the system (for the magnetic field b). It maystill be possible to use the proposed discretisation in space, but the73time discretisation and keeping the consistency of the previously madeassumptions (e.g., non-relativistic motion) may introduce additionalchallenges.6. Other non-linear solvers:The non-linear solver currently used is simply the linearised Oseenform (i.e., a linearly convergent Picard iteration) for the MHD system.A possible next step includes considering other non-linear solvers thathave faster convergence properties. One such scheme is the Newtoniteration, which converges quadratically for a sufficiently accurate ini-tial guess. One challenge here is that the convergence is local, and itis necessary to use various means to compute initial iterates that areclose enough to the solution. A possible approach is to start with afew Picard iterations then use the Newton iterations.7. Different mixed finite element discretisations:We used Taylor-Hood elements [59] for the fluid variables. This choiceof elements is a very common option when considering mixed finiteelement discretisations of fluid flow equations. However, when usingcontinuous pressure elements additional errors may arise due to poormass conservation on the discrete level. A possible example of thisphenomenon has been studied in [39] in the context of colliding flows.One option to overcome this difficulty is to use exactly divergence-free elements such as those proposed in [11, 25]. The elements thereare based on using divergence-conforming elements for the velocitiesand on discontinuous pressure elements. A discontinuous Galerkin74approach is employed to ensure H1-continuity of the velocity fields.The preconditioning approaches in this thesis are based [35], whereexactly divergence-free elements have been used for the fluid unknownsto approximate MHD problem (1.11). One major challenge on thelinear algebra front is the development of a scalable iterative solver forsuch elements. A large-scale study of exactly divergence-free elementsfor the fluid unknowns, and more generally the investigation of othermixed finite element approximations with discontinuous pressures forMHD models is another area of possible future work.75Bibliography[1] M. Alnæs. UFL: a finite element form language. In A. Logg, K. A.Mardal, and G. N. Wells, editors, Automated Solution of DifferentialEquations by the Finite Element Method, Volume 84 of Lecture Notesin Computational Science and Engineering, chapter 17. Springer, 2012.[2] M. Alnæs, A. Logg, K. A. Mardal, O. Skavhaug, and H. P. Langtangen.Unified framework for finite element assembly. International Journalof Computational Science and Engineering, 4(4):231–244, 2009.[3] M. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells. Uni-fied form language: A domain-specific language for weak formulationsof partial differential equations. ACM Transactions on MathematicalSoftware, 40(2):9:1–9:37, 2014.[4] M. S. Alnæs, A. Logg, and K. A. Mardal. UFC: a finite element codegeneration interface. In A. Logg, K. A. Mardal, and G. N. Wells, editors,Automated Solution of Differential Equations by the Finite ElementMethod, Volume 84 of Lecture Notes in Computational Science andEngineering, chapter 16. Springer, 2012.[5] P. R. Amestoy, I. S. Duff, and J. Y. L’Excellent. Multifrontal parallel76distributed symmetric and unsymmetric solvers. Computer Methods inApplied Mechanics and Engineering, 184(2):501–520, 2000.[6] P. R. Amestoy, I. S. Duff, J. Y. L’Excellent, and J. Koster. A fullyasynchronous multifrontal solver using distributed dynamic scheduling.SIAM Journal on Matrix Analysis and Applications, 23(1):15–41, 2001.[7] P. R. Amestoy, A. Guermouche, J. Y. L’Excellent, and S. Pralet. Hybridscheduling for the parallel solution of linear systems. Parallel Comput-ing, 32(2):136–156, 2006.[8] F. Armero and J. C. Simo. Long-term dissipativity of time-steppingalgorithms for an abstract evolution equation with applications to theincompressible MHD and Navier-Stokes equations. Computer Methodsin Applied Mechanics and Engineering, 131(1):41–90, 1996.[9] S. Balay, M. F. Adams, J. Brown, P. Brune, K. Buschelman, V. Ei-jkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes,K. Rupp, B. F. Smith, and H. Zhang. PETSc User’s Manual. TechnicalReport ANL-95/11 - Revision 3.4, Argonne National Laboratory, 2013.[10] S. Balay, M. F. Adams, J. Brown, P. Brune, K. Buschelman,V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C.McInnes, K. Rupp, B. F. Smith, and H. Zhang. PETSc Web page.http://www.mcs.anl.gov/petsc, 2014.[11] B. Cockburn, G. Kanschat, and D. Scho¨tzau. A note on discontinu-ous Galerkin divergence-free solutions of the Navier-Stokes equations.Journal of Scientific Computing, 31(1-2):61–73, 2007.77[12] R. Codina and N. Herna´ndez-Silva. Stabilized finite element approxi-mation of the stationary magneto-hydrodynamics equations. Compu-tational Mechanics, 38(4-5):344–355, 2006.[13] M. Costabel and M. Dauge. Singularities of electromagnetic fields inpolyhedral domains. Archive for Rational Mechanics and Analysis,151(3):221–276, 2000.[14] P. A. Davidson. An Introduction to Magnetohydrodynamics. Cam-bridge Texts in Applied Mathematics. Cambridge University Press,Cambridge, 2001.[15] T. A. Davis. Algorithm 832: UMFPACK v4.3—an unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Soft-ware, 30(2):196–199, June 2004.[16] T. A. Davis. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Soft-ware, 30(2):165–195, June 2004.[17] L. Demkowicz and L. Vardapetyan. Modelling of electromagnetic ab-sorption/scattering problems using hp-adaptive finite elements. Com-puter Methods in Applied Mechanics and Engineering, 152(1):103–124,1998.[18] H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite Elements andFast Iterative Solvers: with Applications in Incompressible Fluid Dy-namics. Oxford University Press, 2005.78[19] R. D. Falgout and U. Yang. hypre: A library of high performancepreconditioners. In Computational Science ICCS 2002, volume 2331of Lecture Notes in Computer Science, pages 632–641. Springer BerlinHeidelberg, 2002.[20] R. Fletcher. Conjugate gradient methods for indefinite systems. InG.Alistair Watson, editor, Numerical Analysis, volume 506 of LectureNotes in Mathematics, pages 73–89. Springer Berlin Heidelberg, 1976.[21] R. W. Freund and N. M. Nachtigal. QMR: a quasi-minimal residualmethod for non-Hermitian linear systems. Numerische Mathematik,60(3):315–339, 1991.[22] J. F. Gerbeau. A stabilized finite element method for the incompressiblemagnetohydrodynamic equations. Numerische Mathematik, 87(1):83–111, 2000.[23] J. F. Gerbeau, C. Le. Bris, and T. Lelie`vre. Mathematical Methods forthe Magnetohydrodynamics of Liquid Metals. Oxford University Press,2006.[24] G. H. Golub and C. Greif. On solving block-structured indefinite lin-ear systems. SIAM Journal on Scientific Computing, 24(6):2076–2092,2003.[25] C. Greif, Li, D. Scho¨tzau, and X. Wei. A mixed finite element methodwith exactly divergence-free velocities for incompressible magnetohy-drodynamics. Computer Methods in Applied Mechanics and Engineer-ing., 199(45-48):2840–2855, 2010.79[26] C. Greif and D. Scho¨tzau. Preconditioners for saddle point linear sys-tems with highly singular (1, 1) blocks. Electronic Transactions onNumerical Analysis, Special Volume on Saddle Point Problems, 22:114–121, 2006.[27] C. Greif and D. Scho¨tzau. Preconditioners for the discretized time-harmonic Maxwell equations in mixed form. Numerical Linear Algebrawith Applications, 14(4):281–297, 2007.[28] M. D. Gunzburger, A. J. Meir, and J. S. Peterson. On the existence,uniqueness, and finite element approximation of solutions of the equa-tions of stationary, incompressible magnetohydrodynamics. Mathemat-ics of Computation, 56(194):523–563, 1991.[29] P. He´non, P. Ramet, and J. Roman. PASTIX: a high-performanceparallel direct solver for sparse symmetric positive definite systems.Parallel Computing, 28(2):301–321, 2002.[30] M. R. Hestenes and E. Stiefel. Methods of conjugate gradients forsolving linear systems. Journal of Research of the National Bureau ofStandards, 49:409–436, 1952.[31] R. Hiptmair and J. Xu. Nodal auxiliary space preconditioning inH(curl) and H(div) spaces. SIAM Journal on Numerical Analysis,45(6):2483–2509, 2007.[32] R. C. Kirby. Algorithm 839: FIAT, a new paradigm for computing finiteelement basis functions. ACM Transactions on Mathematical Software,30(4):502–516, 2004.80[33] R. C. Kirby. FIAT: Numerical construction of finite element basis func-tions,. In A. Logg, K. A. Mardal, and G. N. Wells, editors, AutomatedSolution of Differential Equations by the Finite Element Method, Vol-ume 84 of Lecture Notes in Computational Science and Engineering,chapter 13. Springer, 2012.[34] R. C. Kirby and A. Logg. A compiler for variational forms. ACMTransactions on Mathematical Software, 32(3), 2006.[35] D. Li. Numerical Solution of the Time-Harmonic Maxwell equationsand Incompressible Magnetohydrodynamics Problems. PhD thesis, Uni-versity of British Columbia, 2010.[36] D. Li, C. Greif, and D. Scho¨tzau. Parallel numerical solution of thetime-harmonic Maxwell equations in mixed form. Numerical LinearAlgebra with Applications, 19(3):525–539, 2012.[37] X. S. Li. An overview of SuperLU: Algorithms, implementation,and user interface. ACM Transactions on Mathematical Software,31(3):302–325, 2005.[38] X.S. Li, J. W. Demmel, J. R. Gilbert, L. Grigori, M. Shao, andI. Yamazaki. SuperLU Users’ Guide. Technical Report LBNL-44289,Lawrence Berkeley National Laboratory, 1999. http://crd.lbl.gov/~xiaoye/SuperLU/. Last update: August 2011.[39] A. Linke. Collision in a cross-shaped domain—a steady 2D Navier-Stokes example demonstrating the importance of mass conservation81in CFD. Computer Methods in Applied Mechanics and Engineering,198(41-44):3278–3286, 2009.[40] A. Logg, K. A. Mardal, and G. N. Wells, editors. Automated solution ofdifferential equations by the finite element method, volume 84 of LectureNotes in Computational Science and Engineering. Springer, Heidelberg,2012. The FEniCS book.[41] A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells. FFC: theFEniCS form compiler. In A. Logg, K. A. Mardal, and G. N. Wells,editors, Automated Solution of Differential Equations by the Finite El-ement Method, Volume 84 of Lecture Notes in Computational Scienceand Engineering, chapter 11. Springer, 2012.[42] A. Logg and G. N. Wells. DOLFIN: Automated finite element comput-ing. ACM Transactions on Mathematical Software, 37(2):20:1–20:28,2010.[43] A. Logg, G. N. Wells, and J. Hake. DOLFIN: a C++/Python finiteelement library. In A. Logg, K. A. Mardal, and G. N. Wells, editors,Automated Solution of Differential Equations by the Finite ElementMethod, Volume 84 of Lecture Notes in Computational Science andEngineering, chapter 10. Springer, 2012.[44] P. Monk. Finite Element Methods for Maxwell’s Equations. OxfordUniversity Press, 2003.[45] U. Mu¨ller and L. Bu¨hler. Magnetofluiddynamics in Channels and Con-tainers. Springer, 2001.82[46] M. F. Murphy, G. H. Golub, and A. J. Wathen. A note on precondi-tioning for indefinite linear systems. SIAM Journal on Scientific Com-puting, 21(6):1969–1972, 2000.[47] J. C. Ne´de´lec. Mixed finite elements in R3. Numerische Mathematik,35(3):315–341, 1980.[48] K. B. Ølgaard and G. N. Wells. Optimisations for quadrature represen-tations of finite element tensors through automated code generation.ACM Transactions on Mathematical Software, 37:23, 2010.[49] C. C. Paige and M. A. Saunders. Solution of sparse indefinite systems oflinear equations. SIAM Journal on Numerical Analysis, 12(4):617–629,1975.[50] E. G. Phillips, H. C. Elman, E. C. Cyr, J. N. Shadid, and R. P.Pawlowski. A block preconditioner for an exact penalty formulation forstationary MHD. Technical Report CS-TR-5032, University of Mary-land, Computer Science, 2014.[51] P. H. Roberts. An Introduction to Magnetohydrodynamics. Longmans,London, 1967.[52] Y. Saad. A flexible inner-outer preconditioned GMRES algorithm.SIAM Journal on Scientific Computing, 14(2):461–469, 1993.[53] Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, Philadel-phia, PA, second edition, 2003.83[54] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residualalgorithm for solving nonsymmetric linear systems. SIAM Journal onScientific and Statistical Computing, 7(3):856–869, 1986.[55] D. Scho¨tzau. Mixed finite element methods for stationary incompress-ible magneto–hydrodynamics. Numerische Mathematik, 96(4):771–800,2004.[56] D. J. Silvester and A. J. Wathen. Fast iterative solution of stabilisedStokes systems Part II: Using general block preconditioners. SIAMJournal on Numerical Analysis, 31(5):1352–1367, 1994.[57] T. A. Davis and I. S. Duff. An unsymmetric-pattern multifrontalmethod for sparse LU factorization. SIAM Journal on Matrix Anal-ysis and Applications, 18(1):140–158, 1997.[58] T. A. Davis and I. S. Duff. A combined unifrontal/multifrontal methodfor unsymmetric sparse matrices. ACM Transactions on MathematicalSoftware, 25(1):1–20, 1999.[59] C. Taylor and P. Hood. A numerical solution of the Navier-Stokesequations using the finite element technique. Computers and Fluids,1(1):73–100, 1973.[60] H. A. van der Vorst. Bi-CGSTAB: a fast and smoothly convergingvariant of Bi-CG for the solution of nonsymmetric linear systems. SIAMJournal on Scientific and Statistical Computing, 13(2):631–644, 1992.[61] A. J. Wathen and D. J. Silvester. Fast iterative solution of stabilised84Stokes systems Part I: Using simple diagonal preconditioners. SIAMJournal on Numerical Analysis, 30(3):630–649, 1993.85Appendices86Appendix ACurl operators and crossproductsA.1 2D curlGiven 2D vector fields b(x, y) = (b1, b2), u(x, y) = (u1, u2) and the scalarfunction r(x, y), the curl and cross products are∇× b = ∂b2∂x− ∂b1∂y,∇× r =(∂r∂y,−∂r∂x),u× b =u1b2 − u2b1.Note that taking the curl of a 2D vector field results in a scalar func-tion which is the component in the normal direction to the 2D field (z-component).87A.2 3D curlGiven 3D vector fields b(x, y, z) = (b1, b2, b3) and u(x, y, z) = (u1, u2, u3),the curl and cross products are∇× b =∂b3∂y −∂b2∂z∂b1∂z −∂b3∂x∂b2∂x −∂b1∂y,u× b =u2b3 − u3b2u3b1 − u1b3u1b2 − u2b1.88Appendix BKrylov subspace methodsPDE discretisation usually leads to sparse linear systems. For problems ofsmall dimension, a popular choice to solve such a system are sparse directmethods. However, for systems of large dimension, it is often not possi-ble to use direct methods due to computational memory restrictions, timeconstraints, etc. Therefore an iterative solution method, based on matrix-vector products, is required. One state of the art approach is to considerpreconditioned Krylov subspace methods [53].To review Krylov subspace methods, consider iteratively solving the lin-ear systemAx = b,where A is a real non-singular n×n matrix with vectors x and b of the samedimension. Suppose our initial guess is given by x0. Then the initial residualis given by r0 = b−Ax0. Usually in practice this is taken to be zero. Usingthe initial residual r0, the m-dimensional Krylov subspace is defined as:Km(A; r0) = span{r0, Ar0, A2r0, · · · , Am−1r0}.A Krylov subspace method is based on finding a solution, xm, in the Krylov89subspace:A−1b ≈ xm = x0 + ym,where ym ∈ Km(A; r0) satisfies some optimality criteria. One could charac-terise the three main components to Krylov subspace methods by:1. Computing an orthogonal or bi-orthogonal basis to the Krylov sub-space;2. Defining an optimality criterion within the Krylov subspace;3. Preconditioning approaches.Basis computationsOrthogonal basis:Consider a general real unsymmetric matrix A. Then the process inwhich it is possible to construct an orthogonal basis to the Krylov subspaceKm(A; r0) is known as the Arnoldi process. The Arnoldi process is essen-tially the modified Gram-Schmidt process applied to the Krylov subspace.This process is done iteratively, so at each iteration a new orthogonal vectoris added to the basis. The Arnoldi process is written as:AQm = Qm+1Hm+1,m, (B.1)where Qm+1 is the matrix whose m+1 columns form an orthogonal basis forthe subspace Km(A; r0). Qm has the same leading m columns as in Qm+1and Hm+1,m ∈ R(m+1)×m is an upper Hessenberg matrix. From (B.1) one90can deduce the following form:QTmAQm = Hm,m, (B.2)where Hm,m is the matrix containing the first m rows of Hm+1,m. If A issymmetric then from (B.2) Hm,m must also be symmetric. This leads to theLanzcos method where Hm,m is tridiagonal and hence, we denote it as Tm,min this case. Then the Lanzcos process is as followsAQm = Qm+1Tm+1,m, (B.3)where Tm+1,m has one more extra row than Tm,m. See [53] for more detailsabout both the Arnoldi and Lanczos processes.Bi-Orthogonal basis:For non-symmetric matrices we cannot expect short term recurrences(resulting from the orthogonal basis construction) as well as retain opti-mality in the same way that symmetric solvers do. For example, GMRESminimises the residual but entails long recurrences. The class of biconjugatemethods preserve the short term recurrences property and hence, have lessmemory requirements and are very popular as a result.To introduce the bi-orthogonal basis construction, we define a secondKrylov subspace:Km(AT ; r0) = span{r0, AT r0, (AT )2r0, · · · , (AT )m−1r0}.91Let Vm and Wm be biorthogonal basis matrices with respect to the subspacesKm(A; r0) and Km(AT ; r0), that is W TmVm = I. Then the idea is to performtwo “Lanzcos-type” processes for A and AT . This leads to the procedure:AVm = Vm+1Tˆm+1,m, (B.4a)ATWm = Wm+1Sm+1,m, (B.4b)where Tˆm+1,m and Sm+1,m are tridiagonal matrices of the same structure asTm+1,m as in (B.3). Pre-multiplying (B.4a) by W Tm givesW TmAVm = W TmVm+1Tˆm,m = Tm,m,as W TmVm = I by construction. This is known as (non-symmetric) Lanzcosbiorthogonalisation and leads to such methods as BiCG [20], BiCGstab [60]and QMR [21].Optimality criteriaThere are two popular choices for how to define the optimality criterionfor a Krylov solver. These are minimising the L2-norm within the Krylovsubspace (minimum residual approach) and ensuring that the residual isorthogonal to the Krylov subspace (Galerkin approach).For the three Krylov subspace methods that will be used within this the-sis, namely GMRES [54] for unsymmetric systems, MINRES [49] for sym-metric and Conjugate Gradients [30] (CG) for symmetric positive definitesystems, the optimality criteria are as follows:92• GMRES: min ‖b−Axm‖2, i.e., minimum residual approach• MINRES: min ‖b−Axm‖2, i.e., minimum residual approach• CG: min ‖xm − x‖A, the energy norm of the error which is equiva-lent to orthogonalising with respect to the search space, i.e., GalerkinapproachFor more details about these methods and other Krylov subspace methodswe refer the reader to [53].PreconditioningPreconditioning is an integral part of using any Krylov subspace method.The convergence of CG as well as other methods primarily depend on thespread of the eigenvalues. Preconditioning the linear system Ax = b amountsto solving the system M−1Ax = M−1b, where M is a preconditioner. Thereare certain properties that are required to create a good preconditioner:1. the preconditioner M must approximate A so that the eigenvalues ofM−1A are clustered;2. systems associated with M should be much easier to solve than systemswith A.Loosely speaking, one could classify preconditioners into two types. Thefirst type is an operator-based preconditioner. By this we mean that given acertain problem, usually based on a PDE discretisation, it may be possibleto design certain operators that lead to an optimal or nearly optimal precon-ditioner. These operators are derived from an understanding of the physical93problem and its discretisation. Often these operators are based on spec-tral equivalence: the eigenvalues of M−1A are independent of discretisationparameters. Weak dependency is sometimes sufficient.The second type is more of a “black-box” method, where the precon-ditioner is purely based on the matrix coefficients. The preconditioningstrategies employed in this thesis are mainly based on the first type.94
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Iterative solution of a mixed finite element discretisation...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Iterative solution of a mixed finite element discretisation of an incompressible magnetohydrodynamics… Wathen, Michael 2014
pdf
Page Metadata
Item Metadata
Title | Iterative solution of a mixed finite element discretisation of an incompressible magnetohydrodynamics problem |
Creator |
Wathen, Michael |
Publisher | University of British Columbia |
Date Issued | 2014 |
Description | The aim of this thesis is to develop and numerically test a large scale preconditioned finite element implementation of an incompressible magnetohydrodynamics (MHD) model. To accomplish this, a broad-scope code has been generated using the finite element software package FEniCS and the linear algebra software PETSc. The code is modular, extremely flexible, and allows for implementing and testing different discretisations and linear algebra solvers with relatively modest effort. It can handle two- and three-dimensional problems in excess of 20 million degrees of freedom. Incompressible MHD describes the interaction between an incompressible electrically charged fluid governed by the incompressible Navier-Stokes equations coupled with electromagnetic effects from Maxwell’s equations in mixed form. We introduce a model problem and a mixed finite element discretisation based on using Taylor-Hood elements for the fluid variables and on a mixed N ́ed ́elec pair for the magnetic unknowns. We introduce three iteration strategies to handle the non-linearities present in the model, ranging from Picard iterations to completely decoupled schemes. Adapting and extending ideas introduced in [Dan Li, Numerical Solution of the Time-Harmonic Maxwell Equations and Incompressible Magnetohydrodynamics Problems, Ph.D. Dissertation, The University of British ii Columbia, 2010], we implement a preconditioning approach motivated by the block structure of the underlying linear systems in conjunction with state of the art preconditioners for the mixed Maxwell and Navier-Stokes subproblems. For the Picard iteration scheme we implement an inner-outer preconditioner. The numerical results presented in this thesis demonstrate the efficient performance of our preconditioned solution techniques and show good scalability with respect to the discretisation parameters. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2014-08-26 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivs 2.5 Canada |
DOI | 10.14288/1.0135538 |
URI | http://hdl.handle.net/2429/50202 |
Degree |
Master of Science - MSc |
Program |
Computer Science |
Affiliation |
Science, Faculty of Computer Science, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2014-09 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/2.5/ca/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2014_september_wathen_michael.pdf [ 261.99kB ]
- Metadata
- JSON: 24-1.0135538.json
- JSON-LD: 24-1.0135538-ld.json
- RDF/XML (Pretty): 24-1.0135538-rdf.xml
- RDF/JSON: 24-1.0135538-rdf.json
- Turtle: 24-1.0135538-turtle.txt
- N-Triples: 24-1.0135538-rdf-ntriples.txt
- Original Record: 24-1.0135538-source.json
- Full Text
- 24-1.0135538-fulltext.txt
- Citation
- 24-1.0135538.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0135538/manifest