Scalable implicit incompressible resistive MHD with stabilized FE and fully-coupled Newton–Krylov-AMG☆
Introduction
The resistive magnetohydrodynamics (MHD) model describes the dynamics of charged fluids in the presence of electromagnetic fields. MHD models describe important phenomena in the natural world (e.g., stellar interiors, solar flares, and planetary magnetic field generation) and in technological applications (e.g., magnetically confined fusion energy [e.g. ITER—tokamak], and pulsed fusion reactors such as Z-pinch-type devices) [1]. The mathematical basis for the continuum modeling of these systems is the solution of the governing partial differential equations (PDEs) describing conservation of mass, momentum, and energy, augmented by the low-frequency Maxwell’s equations. This PDE system is non-self adjoint, strongly coupled, highly nonlinear, and characterized by physical phenomena that span a very large range of length- and time-scales. These interacting, nonlinear multiple time-scale physical mechanisms can balance to produce steady-state behavior, nearly balance to evolve a solution on a dynamical time-scale that is long relative to the component time-scales, or can be dominated by a few fast modes. These characteristics make the scalable, robust, accurate, and efficient computational solution of these systems, over relevant dynamical time-scales of interest, extremely challenging.
In the context of MHD simulations, the dominant computational solution strategies for many years have been the use of explicit [2], [3], [4], [5] and partially implicit time integrations methods such as implicit–explicit [6], [7], [8], [9], [10], semi-implicit [11], [12], [13], [14], and operator-splitting [15], [16] techniques. With the exception of fully-explicit strategies, which are limited by numerical stability to follow the fastest normal-mode time-scale, all these temporal integration methods include some implicitness to enhance efficiency by removing one or more sources of numerical stiffness in the problem, either from diffusion or from fast-wave phenomena. While these types of techniques currently form the basis for many production-level resistive MHD simulation tools (see e.g. [14], [9], [17], [16], [18]), a number of applied mathematics, numerical mathematics, and computational science issues remain. Chief among these are the lack of numerical robustness (due to conditional stability constraints stemming from their partially implicit nature), temporal accuracy (which also limits the usable timestep), and the lack of robust, efficient, scalable nonlinear iterative solution methods (which limits the achievable resolution).
From this point of view, fully-implicit formulations, coupled with effective robust nonlinear iterative solution methods, become attractive, as they have the potential to provide stable, higher-order time-integration of these complex multiphysics systems when long dynamical time-scales are of interest. These methods can follow the dynamical time-scales of interest, as opposed to time-scales determined by either numerical stability or by temporal order-of-accuracy reduction. However, as described above, the existence of fast normal modes in MHD make the associated algebraic systems very stiff, and thus difficult to solve iteratively using modern iterative techniques. Despite these challenges, much effort has been invested in the last decade towards the development of more fully-implicit, nonlinearly coupled solution methods for MHD, with the goal of enhancing robustness, accuracy, and efficiency (see e.g. [15], [7], [6], [19], [20], [21], [22], [23], [24], [25], [26], [27]).
In addition to the challenges associated with robustly and efficiently dealing with the myriad of time-scales produced by the resistive MHD system, there are a number of spatial discretization issues that add considerable complexity to the numerical approximation of these systems. These include dealing with the dual saddle-point structure of the velocity–pressure interactions in the Navier–Stokes operators, and the satisfaction/enforcement of the solenoidal involution/constraint on the magnetic induction in conjunction with the evolution of the induction equation for . Briefly, in the context of finite volume (FV) and finite element (FE) methods, approaches to deal with these difficulties include: physics-compatible discretizations that enforce directly key mathematical properties of the continuous problem (see e.g. [28], [29], [30], [31], [32], [22], [33]); methods that transform to potential-based formulations to eliminate one or both saddle-point sub-systems [20], [34], [9], [35], [25], exact and weighted-exact penalty formulations [36], [37], [38], [39], and stabilization methods that consider the Lagrange multiplier form of the solenoidal involution/constraint for the induction equation, and then regularize the double-saddle point structure [40], [41], [42].
This study pursues fully-implicit formulations in the context of a specific FE discretization based on regularization/stabilization methods. These methods introduce a regularization of the magnetic field saddle-point sub-problem to enforce the solenoidal constraint on the magnetic field , and generate a system with four unknowns, the three components of and the Lagrange multiplier, [40], [41], [42]. This system, combined with the Navier–Stokes sub-system, produces a double-saddle point structure for the resulting MHD system and thus requires compatible spaces for both and and that respect the Ladyzhenskaya–Babuska–Brezzi (LLB) condition. The condition for is analogous to the well-known constraint on the velocity and pressure spaces in the Stokes flow system (see e.g. [43], [44]). In this context, references [41], [42] have developed a stabilized FE formulation for resistive MHD in the curl form of the equations with constant properties, and developed weak forms that regularize the double saddle-point structure, and generate streamline upwind Petrov–Galerkin (SUPG) like terms that control oscillation due to convection effects. These authors present a coercivity result for the system that enables development of stabilization parameters that obtain optimal convergence characteristics. It should be noted that an important limitation of these approaches is that they can be shown to be convergent only when the magnetics sub-problem is properly posed in , which applies when the domain is convex or its boundary is sufficiently smooth (see e.g. [37], [45]). Recently, in [45], a stabilized FE method has been proposed that is unconditionally stable and convergent. This stabilized formulation differs from the approaches above in that it defines separate residuals, stabilization terms, and stabilization parameters for the induction and the Lagrange multiplier and requires use of meshes with particular macro-element structures. With these characteristics, convergence is demonstrated on non-convex domains and to solutions with singularities.
The current study complements previous work by developing a robust, efficient, fully-coupled stabilized FE formulation for 3D resistive MHD with solution methods that enable both scalable fully-implicit transient and direct-to-steady-state solutions. This work is an extension of [25] that demonstrated for the first time a scalable solution method for a 2D vector potential formulation for fully-implicit visco-resistive incompressible MHD based on stabilized FE. In that study, the induction equation saddle-point problem is eliminated by use of the 2D vector potential [25]. The current study presents a full 3D formulation that deals with the double saddle-point structure of the incompressible visco-resistive MHD system and demonstrates robust and scalable iterative nonlinear solutions based on fully-coupled algebraic multilevel preconditioned Newton–Krylov methods. It should be noted that the proposed stabilized FE formulation has the additional benefit of enabling the use of equal-order interpolation for all quantities. This simplifies the data structures of the parallel unstructured FE code, and the linear algebra interface for the iterative solution methods that are employed.
The remainder of this paper is organized as follows. Section 2 presents our formulation of the resistive incompressible MHD model. The stabilized FE formulation of the governing 3D resistive MHD equations is developed by a variational multiscale (VMS) method and presented in Section 3. In Section 4, a brief overview of the fully-implicit Newton–Krylov solution method is presented with a discussion of the domain decomposition and multilevel preconditioners. Section 5 presents representative order-of-accuracy verification studies and performance, scaling and simulation results for some illustrative resistive MHD problems. Finally, Section 6 closes with a few conclusions.
Section snippets
Resistive MHD model equations
The system of interest for this study is the 3D resistive iso-thermal MHD equations including dissipative terms for the momentum and magnetic induction equations [1]. This model provides a continuum description of charged fluids in the presence of electromagnetic fields. The system of equations is: Here is the plasma velocity; is the ion mass density; is the stress tensor; is the magnetic induction (here after also
A stabilized finite element formulation for 3D resistive MHD
The continuous PDE problem, defined by the 3D resistive MHD equations in Table 1, is approximated by a stabilized FE formulation. This formulation allows for stable equal-order velocity–pressure interpolation and provides for convection stabilization, as described below.
We employ stabilized FE methods to avoid stability and algorithmic limitations of mixed Galerkin FE formulations. In particular, in a mixed Galerkin FE formulation of the momentum-continuity equations of the Navier–Stokes part
Fully-implicit time integration and direct to steady-state solutions
For multiple-time-scale multiphysics PDE systems such as resistive MHD, fully-implicit methods are an attractive choice that can often provide unconditionally-stable time integration techniques. These methods can be designed with various types of stability properties that allow robust integration of multiple timescale systems without the requirement to resolve the stiff modes of the system (which are not of interest when they do not control the accuracy of time integration [66]). In this study
Results and discussion
In this study, we present two general classes of results. The first in Section 5.1, is for assessment of the accuracy and order-of-accuracy of the stabilized FE formulation described above by a comparison with analytic solutions. The problems in this section include, (1) a steady-state MHD duct flow with strong coupling between the momentum and the Lorentz forces that provides a detailed assessment of the order-of-accuracy of the spatial discretization; and (2) a transient MHD Rayleigh flow
Conclusions
This paper has presented a new scalable fully-implicit 3D resistive MHD formulation based on an unstructured stabilized FE formulation. The formulation was briefly presented as developed by a variational multiscale (VMS) development of the divergence form of the governing balance equations. The formulation employed a scalar Lagrange multiplier that was used to enforce the solenoidal involution for the magnetic field as a constraint. A critical aspect of the formulation is the stabilization of
Acknowledgments
The authors would like to thank Paul Lin for continued help with improvement of the usability and scaling of the Trilinos suite of solution methods software on large-scale HPC platforms, and David Sondak for suggestions and comments on an early version of the manuscript. The parallel scaling studies were carried out on the DOE/Oak Ridge National Laboratory Jaguar CrayXK7 platform and performed under computing resource allocation from the Joule Metric Effort.
References (101)
- et al.
A simple finite difference scheme for multidimensional magnetohydrodynamic equations
J. Comput. Phys.
(1998) Divergence-free adaptive mesh refinement for magnetohydrodynamics
J. Comput. Phys.
(2001)The constraint in shock-capturing magnetohydrodynamics codes
J. Comput. Phys.
(2000)- et al.
An implicit algorithm for compressible three-dimensional magnetohydrodynamic calculations
J. Comput. Phys.
(1985) - et al.
Semi-implicit method for three-dimensional compressible magnetohydrodynamic simulation
J. Comput. Phys.
(1985) - et al.
Semi-implicit magnetohydrodynamic calculations
J. Comput. Phys.
(1987) - et al.
Accurate semi-implicit treatment of the Hall effect in magnetohydrodynamic computations
J. Comput. Phys.
(1989) - et al.
Nonlinear magnetohydrodynamics simulation using high-order finite elements
J. Comput. Phys.
(2004) - et al.
Implicit, nonlinear reduced resistive MHD nonlinear solver
J. Comput. Phys.
(2002) - et al.
An implicit nonlinear reduced resistive MHD solver
J. Comput. Phys.
(2002)
Approximate riemann solver for the two-fluid plasma model
J. Comput. Phys.
A non-staggered, conservative, , finite-volume scheme for 3D implicit extended magnetohydrodynamics in curvilinear geometries
Comput. Phys. Comm.
Towards a scalable fully-implicit fully-coupled resistive MHD formulation with stabilized FE methods
J. Comput. Phys.
Advanced physics calculations using a multi-fluid plasma model
Comput. Phys. Comm.
Adjoint operators for the natural discretizations of the divergence, gradient and curl on logically rectangular grids
Appl. Numer. Math.
Hall MHD effects in the 2-d Kelvin-Helmholtz/tearing instability
Phys. Lett. A
An adaptive finite element method for magnetohydrodynamics
J. Comput. Phys.
Approximation of the thermally coupled mhd problem using a stabilized finite element method
J. Comput. Phys.
On an unconditionally convergent stabilized finite element approximation of resistive magnetohydrodynamics
J. Comput. Phys.
Hyperbolic divergence cleaning for the mhd equations
J. Comput. Phys.
Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations
Comput. Methods Appl. Mech. Engrg.
A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuska-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations
Comput. Methods Appl. Mech. Engrg.
Multiscale phenomena: Green’s functions, the dirichlet-to-neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods
Comput. Methods Appl. Mech. Engrg.
The variational multiscale method: A paradigm for computational mechanics
Comput. Methods Appl. Mech. Engrg.
Stabilized finite element approximation of transient incompressible flows using orthogonal subscales
Comput. Methods Appl. Mech. Engrg.
A new class of finite element variational multiscale turbulence models for incompressible magnetohydrodynamics
J. Comput. Phys.
A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimentional advective-diffusive systems
Comput. Methods Appl. Mech. Engrg.
The effect of nonzero on the numerical solution of the magnetohydrodynamic equations
J. Comput. Phys.
A method for incorporating Gauss law into electromagnetic pic codes
J. Comput. Phys.
An inexact Newton method for fully-coupled solution of the Navier–Stokes equations with heat and mass transport
J. Comput. Phys.
Stabilized FE computational analysis of nonlinear steady state transport/reaction systems
Comput. Methods Appl. Mech. Engrg.
Performance of fully-coupled domain decomposition preconditioners for finite element transport/reaction simulations
J. Comput. Phys.
Hall MHD effects in the 2-D Kelvin-Helmholtz/tearing instability
Phys. Lett.: A
Principles of Magnetohydrodynamics with Applications to Laboratory and Astrophysical Plasmas
A divergence-free upwind code for multi-dimensional magnetohydrodynamics flows
Astrophys. J.
Implicit and semi-implicit schemes in the versatile advection code: numerical tests
Astronom. Astrophys.
Implicit and semi-implicit schemes: algorithms
Internat. J. Numer. Methods Fluids
Nonlinear simulation studies of tokamaks and STS
Nucl. fusion
Implicit solution of the four-field extended-magnetohydrodynamic equations using high-order high-continuity finite elements
Phys. Plasmas
IRMHD: an implicit radiative and magnetohydrodynamical solver for self-gravitating systems
Mon. Not. R. Astron. Soc.
Three-dimensional hydra simulations of national ignition facility targets
Phys. Plasmas
Coalescence of magnetic islands in the low-resistivity, Hall-MHD regime
Phys. Rev. Lett.
An optimal, parallel, fully implicit newton-krylov solver for three-dimensional visco-resistive magnetohydrodynamics
Phys. Plasmas
Mixed finite elements in
Numer. Math.
Simulation of magnetohydrodynamic flows: a constrained transport method
Astrophys. J.
Matching algorithms with physics: exact sequences of finite element spaces
Cited by (65)
On a fully-implicit VMS-stabilized FE formulation for low Mach number compressible resistive MHD with application to MCF
2023, Computer Methods in Applied Mechanics and EngineeringStructure-preserving and helicity-conserving finite element approximations and preconditioning for the Hall MHD equations
2023, Journal of Computational PhysicsA decoupled, unconditionally energy-stable and structure-preserving finite element scheme for the incompressible MHD equations with magnetic-electric formulation
2023, Computers and Mathematics with ApplicationsStabilized bi-cubic Hermite Bézier finite element method with application to gas-plasma interactions occurring during massive material injection in tokamaks
2023, Computers and Mathematics with ApplicationsA multilevel block preconditioner for the HDG trace system applied to incompressible resistive MHD
2023, Computer Methods in Applied Mechanics and EngineeringCitation Excerpt :The potential benefit of implicit methods however requires robust, efficient, and scalable nonlinear and linear iterative solvers/preconditioners which can handle the highly stiff algebraic systems generated from the incompressible MHD equations. Over the past few years there has been significant progress in the area of fully implicit robust and scalable MHD simulations as evidenced in [2–8] among others. The vast majority of current methods employ low-order stabilized FEM, mixed FEM or finite volume methods for the spatial discretization and hence require larger numbers of unknowns for high accuracy.
A conservative finite element solver for the induction equation of resistive MHD: Vector potential method and constraint preconditioning
2022, Journal of Computational Physics
- ☆
This work was partially supported by DOE NNSA ASC Algorithms effort, the DOE Office of Science AMR program at Sandia National Laboratory under contract DE-AC04-94AL85000, and at Los Alamos National Laboratory under contract DE-AC52-06NA25396.