An a posteriori error analysis for the equations of stationary incompressible magnetohydrodynamics
AAn a posteriori error analysis for the equations ofstationary incompressible magnetohydrodynamics
J. H. Chaudhry a , A. E. Rappaport b,a, ∗ , J. N. Shadid b,a a Department of Mathematics and Statistics, University of New Mexico, Albuquerque NM,87123 b Center for Computing Research, Sandia National Laboratories, Albuquerque NM, 87185
Abstract
Magnetohydrodynamics (MHD) is a continuum level model for conducting fluidssubject to external magnetic fields, e.g. plasmas and liquid metals. The efficientand robust solution of the MHD system poses many challenges due to it’s non-linear, non self-adjoint, and highly coupled nature. In this paper, we developa robust and accurate a posteriori error estimate for the numerical solution ofthe MHD equations based on the exact penalty method. The error estimatealso isolates particular contributions of error in a quantity of interest (QoI) toinform discretization choices to arrive at accurate solutions. The tools requiredfor these estimates involve duality arguments and computable residuals.
1. Introduction
The magnetohydrodynamics (MHD) equations provide a continuum modelfor conducting fluids subject to magnetic fields and are often used to modelimportant applications e.g. highly collisional plasmas. In this context, MHDcalculations aid physicists in understanding both thermonuclear fusion and as-trophysical plasmas as well as understanding the behavior of liquid metals [1, 2].The governing equations of MHD couple Navier-Stokes equations for fluid dy-namics with a reduced set of Maxwell’s equations for low frequency electromag-netic phenomenon. Structurally, the equations of MHD form highly coupled,nonlinear, non self-adjoint partial differential equations (PDEs). Analytical so-lutions to the MHD system cannot be obtained for practical configurations;instead numerical solutions are sought. Solution methods for the incompress-ible MHD system include stabilized formulations, fully implicit block precon-ditioning, First Order System Least Squares and operator splitting methods[3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. In this article we restrict ∗ Corresponding author
Email addresses: [email protected] (J. H. Chaudhry), [email protected],[email protected] (A. E. Rappaport), [email protected] (J. N. Shadid)
Preprint submitted to Elsevier April 27, 2020 a r X i v : . [ m a t h . NA ] A p r urselves to the stationary MHD equations based on the exact penalty formu-lation [3].The numerical solution of complex equations like the MHD equations of-ten have a significant discretization error. This error must be quantified forthe reliable use of MHD equations in numerous science and engineering fields.Accurate error estimation is a key component of predictive computational sci-ence and uncertainty quantification [18, 19, 20]. Moreover, the error dependson a complex interaction between many contributions. Thus, the availabilityof accurate error estimation also offers the potential of optimizing the choiceof discretization parameters in order to achieve desired accuracy in an efficientfashion. In this work we leverage adjoint based a posteriori error estimates fora QoI related to to the solution of the MHD equations. These estimates providea concrete error analysis of different contributions of error, as well as informsolver and discretization strategies.In many scientific and engineering applications, the goal of running a sim-ulation is to compute certain a quantity of interest (QoI) of the solution, forexample the drag over a plane wing in a compressible CFD context. Adjointbased analysis for quantifying the error in a numerically computed QoI has foundsuccess for a wide variety of numerical methods and discretizations ranging fromfinite element, finite difference, finite volume, time integration, operator split-ting techniques and uncertainty quantification [21, 22, 23, 24, 25, 26, 27, 28,29, 30, 31, 32, 33, 34, 35, 20]. Adjoint based a posteriori error analysis usesvariational analysis and duality to relate errors to computable residuals. Inparticular, one solves an adjoint problem whose solution provides the residualweighting to produce the error in the QoI. The technique also naturally allowsto identify and isolate different components of error arising from different as-pects of discretization and solution methods, by analyzing different componentsof the weighted residual separately.This paper is organized as follows. In §
2, we review the equations of in-compressible resistive MHD, present the exact penalty weak form and the finiteelement method to numerically solve the problem. In § a posteriori error analysis for an abstract problemrepresentative of the exact penalty weak form. We apply these results to theMHD equations in § a posteriori error estimate. In § §
2. Exact penalty formulation and discretization for incompressibleMHD
In this section we describe the nondimensionalized equations of incompress-ible stationary MHD, a stabilized weak form of the MHD system and a finiteelement method for its solution. 2 .1. The MHD equations
We consider a Lipschitz domain Ω ⊂ R d , d =2 or 3, with boundary ∂ Ω. Theequations for stationary incompressible MHD in Ω are given by − R f ∆ u + u · ∇ u + ∇ p − κ ( ∇ × b ) × b = f , (2.1a) ∇ · u = 0 , (2.1b) κR m ∇ × ( ∇ × b ) − κ ∇ × ( u × b ) = , (2.1c) ∇ · b = 0 , (2.1d)where the unknowns are the velocity u , the magnetic field b , and the pres-sure p . The nondimensionalized parameters are the fluid Reynolds number R f ,Magnetic Reynolds number R m , and interaction parameter κ = H a ( R f R m ),where H a is the Hartmann number. The source term f is viewed as data to theproblem. For x ∈ Ω we have u ( x ) ∈ R d , b ( x ) ∈ R d , p ( x ) ∈ R and f ( x ) ∈ R d .We supplement the system (2.1) with boundary conditions, u = g , on ∂ Ω , (2.2a) b × n = q × n , on ∂ Ω . (2.2b)Referring to (2.1), we observe there are 2 d + 2 and only 2 d + 1 unknowns [5].Effectively enforcing the solenoidal constraint (2.1d) (an involution of the tran-sient MHD system) is an open area of research. Techniques include compatiblediscretizations [36, 37], vector potential [38, 4] and divergence cleaning [39, 40]as well as the exact penalty method [3, 41, 5]. In this article, we consider theexact penalty method which we further describe in § We make use of the standard spaces L (Ω) and H m (Ω) as well as theirvector counterparts L (Ω) and H m (Ω). The L (Ω) (or L (Ω)) inner product isdenoted by ( · , · ) and the norm is denoted by (cid:107) · (cid:107) , while the H (Ω) (or H (Ω))norm is denoted by (cid:107) · (cid:107) . These details of these function spaces are givenin Appendix A. For b ∈ H (Ω), we define ∇ b := (cid:2) ∇ b , . . . , ∇ b d (cid:3) as a matrixwhose columns are the gradients of the components of b . The relevant subspacesof H (Ω) needed to satisfy the boundary conditions (in the sense of the traceoperator) are, H (Ω) := { w ∈ H : w | ∂ Ω ≡ } , (2.3) H τ (Ω) := { w ∈ H : ( w × n ) | ∂ Ω ≡ } . (2.4)Finally, we define the product space, M (Ω) := H (Ω) × H τ (Ω) × L (Ω) . (2.5)3e also remark that for d = 2, we use the natural inclusion of R (cid:44) → R , (cid:2) v , v (cid:3) T (cid:55)→ (cid:2) v , v , (cid:3) T to define the operators ∇× and × . Thus for v , w ∈ H ,we have that ∇ × v = (cid:18) ∂v y ∂x − ∂v x ∂y (cid:19) ˆ k , v × w = ( v x w y − v y w x ) ˆ k . In this section we present the weak form of the stationary incompressibleMHD system based on the exact penalty formulation. The exact penalty methodrequires that the domain Ω is bounded, convex and polyhedral. This ensuresthat H ( curl , Ω) ∩ H (div , Ω) is continuously embedded in H (Ω) [37]. Wealso assume homogeneous Dirichlet conditions, g = q = . Non-homogeneousboundary conditions can be dealt with through standard lifting arguments asdiscussed in § U = ( u , b , p ) ∈ M (Ω) such that N EP ( U, V ) = ( f , v ) , ∀ V ∈ M (Ω) , (2.6)where the nonlinear form N EP is defined for all V = ( v , c , q ) ∈ M (Ω) by N EP ( U, V ) : = 1 R ( ∇ u , ∇ v ) + ( C ( u ) , v ) − ( p, ∇ · v ) + ( q, ∇ · u )+ κ ( Y ( b ) , v ) − κ ( Z ( u , b ) , c )+ κR m ( ∇ × b , ∇ × c ) + κR m ( ∇ · b , ∇ · c ) , (2.7)and the nonlinear operators are defined by C ( u ) := u · ∇ u , (2.8a) Y ( b ) := ( ∇ × b ) × b , (2.8b) Z ( u , b ) := ∇ × ( u × b ) . (2.8c)All except the last term in the weak form arise from multiplying (2.1a)-(2.1c)by test functions and performing integration by parts. The last term, κR m ( ∇ · b , ∇ · c ), effectively enforces the solenoidal involution (2.1d) since, assuming theaforementioned restrictions on the domain, there exists a function (see [3, 42]) b ∈ H (Ω) such that ∇ · ∇ b = ∇ · b , and ∇ b ∈ H τ (Ω) . (2.9)Thus, we choose V = ( , ∇ b ,
0) in (2.7) and use (B.1d) so that (2.6) reducesto ( ∇ · b , ∇ · ∇ b ) = ( ∇ · b , ∇ · b ) = 0 , (2.10)and hence (2.1d) is satisfied almost everywhere in Ω.4 .4. Finite element method We introduce the standard continuous Lagrange finite element spaces. Let T h be a simplicial decomposition of Ω, where h denotes the maximum diameterof the elements of T h . The standard Lagrange space finite element space of order q is then P qh := (cid:8) v ∈ C (Ω) : ∀ K ∈ T h , v | K ∈ P q ( K ) (cid:9) , (2.11)where P q ( K ) is the space of polynomials of degree at most q defined on theelement K . Additionally, our finite element space satisfies the Ladyzhenskaya-Babuˇska-Brezzi condition stability condition [43] for the velocity pressure pair,e.g. M h (Ω) = P h (Ω) × P h (Ω) × P h (Ω). Then the discrete problem to find anapproximate solution U h = ( u h , b h , p h ) ∈ M h (Ω) to (2.7) is, N EP ( U h , V h ) = ( f , v h ) ∀ V h ∈ M h (Ω) . (2.12)Note there is no restriction on the finite element space for b h , which is anadvantage of this method. The well-posedness of the continuous and discreteproblems are are proven in [3]. The goal of a numerical simulation is often to compute some functional ofthe solution, that is, the QoI. In particular, QoIs considered in this article havethe generic form, QoI = (cid:90) Ω Ψ · U d x = (Ψ , U ) (2.13)where Ψ ∈ L (Ω) × L (Ω) × L (Ω) ≡ [ L (Ω)] d +1 . For example in 2D, tocompute the average of the y component of velocity u y over a region Ω c ⊂ Ω,set Ψ = | Ω c | (0 , Ω c , , , T , where S denotes the characteristic function over aset S . In the examples presented later, the QoIs physically represent quantitiesrepresentative of the average flow rate, or the average induced magnetic field.We seek to compute error estimates in the QoI using duality arguments aspresented in the following subsection.
3. Abstract a posteriori error analysis
In this section we consider an abstract variational setting for a posteriori analysis based on the ideas from [21, 31, 22, 23, 24]. Given a Banach space V ,we denote its dual space by V (cid:48) and represent the duality pairing by (cid:104)· , ·(cid:105) . Weconsider generic QoI as bounded linear functionals of the form, Q ( w ) = (cid:104) ψ, w (cid:105) , where ψ ∈ V (cid:48) and w ∈ V . For example, in (2.13), (cid:104) ψ, u (cid:105) = (Ψ , U ), the dualitypairing is represented by the L inner product. We want to evaluate Q ( u ) where u is the solution to the variational problem: find u in V such that N ( u, v ) = (cid:104) f, v (cid:105) , ∀ v ∈ V , (3.1)5nd N : V × V → R is linear in the second argument but may be nonlinearin the first argument. Throughout this section u refers to the true solution to(3.1). An example of such a variational problem is the exact penalty problemas described in § u h ∈ V h , in some finitedimensional subspace V h ⊂ V , we define the error as e = u − u h . The aim ofthe a posteriori analysis is to compute the error in the QoI, Q ( u ) − Q ( u h ) = (cid:104) ψ, u (cid:105) − (cid:104) ψ, u h (cid:105) = (cid:104) ψ, e (cid:105) . For nonlinear forms, the choice of an adjoint formis not straightforward. However, a common choice useful for various kinds ofanalysis is based on linearization [44, 45, 46, 47, 32, 30]. This choice enables thedefinition of a bilinear form N ∗ ( v, w ) which satisfies the useful property, N ∗ ( v, e ) = N ( u, v ) − N ( u h , v ) , (3.2)for all v ∈ V .Now let V = (cid:81) ni =1 V i be a product space of Banach spaces V i = V i (Ω), withΩ ⊂ R d . The problem (3.1) with N : V × V → R is now more specifically givenby N ( u, v ) = m (cid:88) i =1 (cid:104) N i ( u ) , v (cid:96) i (cid:105) + a ( u, v ) , (3.3)where a ( u, v ) is a bilinear form and (cid:96) i ∈ { , . . . , n } and N i : V → V (cid:48) (cid:96) i . For asolution/approximation pair ( u/u h ), define the matrix J , by J ij = (cid:90) ∂N i ∂u j ( su + (1 − s ) u h ) d s, (3.4)where ∂N i ∂u j ( · ) denotes the partial derivative of N i with respect to the argument u j . Define the linearized operator ¯ N i for w ∈ V by¯ N i w = (cid:90) ∂N i ∂u ( su + (1 − s ) u h ) d s · w = n (cid:88) j =1 (cid:90) ∂N i ∂u j ( su + (1 − s ) u h ) d s w j = n (cid:88) j =1 J ij w j . Now since each ¯ N i is linear, we may define the bilinear forms, ν i ( u, v ) = (cid:104) ¯ N i u, v (cid:96) i (cid:105) = (cid:42) n (cid:88) j =1 J ij u j , v (cid:96) i (cid:43) = n (cid:88) j =1 (cid:10) J ij u j , v (cid:96) i (cid:11) . Define ν ∗ i ( v, w ) = ν i ( w, v ), adjoint operators J ∗ ij to J ij satisfying (cid:104)J ij w, v (cid:105) = (cid:104) w, J ∗ ij v (cid:105) (3.5)for w ∈ V j and v ∈ V (cid:48) l i and a ∗ ( w, v ) := a ( v, w ). With these definitions, theadjoint bilinear weak form is, N ∗ ( φ, w ) = m (cid:88) i =1 ν ∗ i ( φ, v ) + a ∗ ( φ, v ) = m (cid:88) i =1 n (cid:88) j =1 (cid:104) v j , J ∗ ij φ (cid:96) i (cid:105) + a ∗ ( φ, v ) . (3.6)6hen if φ solves the dual problem, N ∗ ( φ, v ) = (cid:104) ψ, v (cid:105) , ∀ v ∈ V , (3.7)we then have the following abstract error representation. Theorem 1.
The error in a QoI represented by Q ( u ) = (cid:104) ψ, u (cid:105) is computable as (cid:104) ψ, e (cid:105) = (cid:104) f, v (cid:105) − N ( u h , φ ) .Proof. Unpacking the definitions, (cid:104) ψ, e (cid:105) = N ∗ ( φ, e ) = m (cid:88) i =1 n (cid:88) j =1 (cid:104) e j , J ∗ ij φ (cid:96) i (cid:105) + a ∗ ( φ, e )= m (cid:88) i =1 n (cid:88) j =1 (cid:104)J ij e j , φ (cid:96) i (cid:105) + a ( e, φ ) = m (cid:88) i =1 (cid:104) N i e, φ (cid:96) i (cid:105) + a ( e, φ )= m (cid:88) i =1 (cid:104) N i ( u ) − N i ( u h ) , φ (cid:96) i (cid:105) + a ( u, φ ) − a ( u h , φ )= m (cid:88) i =1 (cid:104) N i ( u ) , φ (cid:96) i (cid:105) + a ( u, φ ) − m (cid:88) i =1 (cid:104) N i ( u h ) , φ (cid:96) i (cid:105) − a ( u h , φ )= N ( u, φ ) − N ( u h , φ ) = ( f, φ ) − N ( u h , φ ) . The main takeaway of this theorem is that computing the adjoint to a non-linear form is reduced to computing the adjoint for the averaged entries, J ij . A posteriori error estimate for the MHD equations
The analysis in § L (Ω)] d +1 inner product ( · , · ). Thelinear and nonlinear terms in the exact penalty weak form (2.6) are mapped tomatch (3.3). The mapping between the abstract formulation and MHD equationis shown in Table 1.For the exact penalty weak form, we have that N EP ( U, V ) = (cid:88) i =1 ( N EP,i ( U ) , V (cid:96) i ) + a EP ( U, V ) , (4.1)where ( N EP, ( U ) , V ) = ( Z ( u , b ) , v ) , ( N EP, ( U ) , V ) = ( Y ( b ) , c ) , ( N EP, ( U ) , V ) = ( C ( u ) , v ) , (4.2)7 bstract MHD (cid:104) , (cid:105) ( , ) m N N EP u Uv VN i N EP,i (a)
Abstract MHD (cid:104) f, v (cid:105) ( f , v ) u U ≡ u u U ≡ b u U ≡ pv V ≡ v v V ≡ c (b) Abstract MHD v V ≡ q J ∗ Z ∗ u J ∗ Z ∗ b J ∗ Y ∗ J ∗ C ∗ a a EP (c)Table 1: Mapping between the abstract framework in § § N EP is given in (4.1), N EP,i in (4.2), a EP in (4.3) and Z ∗ u , Z ∗ b , Y ∗ , C ∗ are given in (4.4). Z , Y , C are in turned defined in (2.8), and a EP ( U, V ) = 1 R ( ∇ u , ∇ v ) − ( p, ∇ · v ) + ( q, ∇ · u )+ κR m ( ∇ × b , ∇ × c ) + κR m ( ∇ · b , ∇ · c ) . (4.3)The entries J ∗ V = Z ∗ u , J ∗ V = Z ∗ b , J ∗ V = Y ∗ v and J ∗ V = C ∗ v are, Z ∗ u c = ( u + u h ) × ( ∇ × c ) , Z ∗ b c = − ( b + b h ) × ( ∇ × c ) , Y ∗ v = (cid:0) − ( ∇ × ( b + b h ) × v ) + ∇ × (( b + b h ) × v ) (cid:1) , C ∗ v = (cid:0) v · (cid:0) ( ∇ u + ∇ u h ) T (cid:1) − (( u + u h ) · ∇ v ) − (( ∇ · ( u + u h ) v ) (cid:1) , (4.4)while the remaining J ∗ ij entries are zero. The details of the derivation are givenin § We are now prepared pose a weak adjoint problem corresponding to exactpenalty primal problem (2.6). Based on (4.1), (4.4) and (3.7), The weak dualproblem is therefore be stated as: find Φ = ( φ , β , π ) ∈ M (Ω) such that N ∗ EP (Φ , V ) = (Ψ , V ) , ∀ V ∈ M (Ω) (4.5)with N ∗ EP (Φ , V ) = 1 R f ( ∇ φ , ∇ v ) + (cid:16) C ∗ φ , v (cid:17) − ( ∇ · v , π ) − ( ∇ · φ , q )+ κR m ( ∇ × β , ∇ × c ) + κR m ( ∇ · β , ∇ · c ) − κ (cid:16) Y ∗ φ , c (cid:17) − κ (cid:16) Z ∗ u β , v (cid:17) − κ (cid:16) Z ∗ b β , c (cid:17) . (4.6)The forms of the linear operators C ∗ , Y ∗ , Z ∗ u and Z ∗ b are given in (4.4). Wediscuss the well-posedness of the adjoint weak form in § .2. Error representation In order to discuss error representation we need to the following definitions
Definition 1.
Define the monolithic error by E = (cid:2) e u , e b , e p (cid:3) T with componenterrors e u = u − u h , e b = b − b h , e p = p − p h . (4.7)We have the following error representation. Theorem 2 (Error representation for exact penalty) . For a given QoI repre-sented by
Ψ = (cid:2) ψ u , ψ b , ψ p (cid:3) T , the error in the numerical approximation of theQoI satisfies (Ψ , E ) = ( f , φ ) − (cid:20) R f ( ∇ u h , ∇ φ ) + ( u h · ∇ u h , φ ) − ( p h , ∇ · φ ) + κ (( ∇ × b h ) × b h , φ ) + ( ∇ · u h , π )+ κR m ( ∇ × b h , ∇ × β ) + κ ( ∇ × ( u h × b h ) , β )+ κR m ( ∇ · b h , ∇ · β ) (cid:21) . Proof.
By Theorem 1,(Ψ , E ) = N ∗ EP (Φ , E ) = N EP ( U, Φ) − N EP ( U h , Φ) = ( f , φ ) − N EP ( U h , Φ) . The analysis above easily extends to the case of non-homogeneous boundaryconditions, i.e. when g or q × n are not identically zero. First assume that thenumerical solution U h on the satisfies the non-homogeneous conditions exactly.That is, u = u h = g and b × n = b h × n = q × n on ∂ Ω. Then, although neitherthe true solution U nor the numerical solution U h belong to M (Ω), yet the error,E = U − U h , satisfies homogeneous boundary conditions and hence belongs to M (Ω). Thus, the error analysis in the previous section applies directly in thiscase.Now, if U h belongs to M h (Ω), then in general U h does not satisfy the non-homogeneous boundary conditions exactly. Hence we consider the splitting ofthe numerical solutions as, U h = U h + U d , (4.8)where U h ∈ M h (Ω) solves, N EP ( U h , V h ) = N EP ( U h + U d , V h ) = ( F, V h ) , ∀ V h ∈ M h (Ω) , (4.9)and U d is a known function that satisfies the non-homogeneous boundary con-ditions accurately. That is, the unknown is now U h and the numerical solution U h is formed through the sum in (4.8). In this article the function U d is approx-imated through a finite element space of much higher dimension than M h (Ω)to capture the boundary conditions accurately.9 .4. Error Estimate and Contributions The error representation in Theorem 2 requires the exact solution Φ =( φ , β , π ) ∈ M (Ω). Moreover, the adjoint weak form (4.6) is linearized aroundthe true solution U and the approximate solution U h . In practice, the adjointsolution itself must be approximated in a finite element space W h ⊂ M (Ω) andis linearized only around the numerical solution. These approximations lead toan error estimate from the error representation 2. Let this estimate be denotedby η . That is, η ≈ (Ψ , E ) where, η = E mom + E con + E M , (4.10)where, E mom = ( f , φ h ) − (cid:18) R f ( ∇ u h , ∇ φ h ) + ( u h · ∇ u h , φ h ) − ( p h , ∇ · φ h )+ κ (( ∇ × b h ) × b h , φ h ) (cid:19) ,E con = − ( ∇ · u h , π h ) ,E M = − κR m ( ∇ × b h , ∇ × β h ) + κ ( ∇ × ( u h × b h ) , β h ) − κR m ( ∇ · b h , ∇ · β h ) . (4.11)Here E mom , E con and E M represent the momentum error contribution, thecontinuity error contribution and the magnetic error contribution respectively.To obtain an accurate error estimate we choose W h to be of much higherdimension than M h (Ω) as is standard in adjoint based a posteriori error esti-mation [25, 21, 31, 48, 35, 25, 49, 50, 34]. Moreover, the inaccuracy caused bysubstituting the numerical solution in place of true solution in the adjoint weakform is shown to decrease in the limit of refined discretization [25].
5. Numerical results
In this section we present numerical results to verify the accuracy of theerror estimate (4.10) and the and utility of the error contributions in (4.11).The effectivity ratio, denoted Eff., characterizes how well the error estimateapproximates the true error,Eff. = Error estimateTrue error = η (Ψ , E ) . (5.1)The closer the effectivity is to 1, the better the error estimate provided by ourmethod.We present two numerical examples here, the Hartmann problem § § Dolfin in the
FEniCS suite [51, 52, 53].10or all experiments, we chose different polynomial orders of Lagrange spacesfor the product space M h (Ω) and ensure that the adjoint space W h has a higherpolynomial degree. The computational domain for all problems is chosen to bea unit length square, Ω := [ − , ] ⊂ R . The mesh is a simplicial uniformmesh with the total number of elements 2D Elem. = n e × n e . Our first results concern the so-called Hartmann problem [2]. This problemmodels the flow of a conducting fluid in a channel. In this case we take considera square channel, as described in the beginning of the section. This problemadmits an analytic solution [5], u = ( u x , b = ( B x , , p where u x ( y ) = G R f (cosh( H a / − cosh( H a y ))2 H a sinh( H a / , (5.2a) B x ( y ) = G (sinh( H a y ) − H a / y )2 κ sinh( H a / , (5.2b) p ( x ) = − Gx − κB x / , (5.2c)and G = − dpdx is an arbitrary pressure drop that we choose to normalize maxi-mum velocity u x ( y ) to 1. The values of the nondimensionalized constants are chosen as follows: R f =16 , R m = 16 , κ = 1 which produce a Hartmann number of H a = 16. The QoI ischosen as the average velocity across the flow over a slice. To this end, defineΩ c := (cid:2) − , (cid:3) × (cid:2) − , (cid:3) (5.3)and consequently Ω c the characteristic function on Ω c . The monolithic quantityof interest Ψ as in Theorem 2 is chosen to be Ψ = ( Ω c , , , , T . The QoIthus reduces to (Ψ , U ) = ( Ω c , u x ) . (5.4)This has a physical interpretation of the capturing the flow rate across this sliceof the channel, Ω c . Error contributions and effectivites using different order polynomial spacesare presented in Table 2, Table 3, Table 4, and Table 5. Effectivity ratio intables Table 2 and Table 3 is quite close to 1 indicating the accuracy of the errorestimate. The error estimate in Table 4 is not as accurate due to linearizationerror incurred by replacing the true solution by the approximate solution inthe definition of the adjoint as discussed in the introduction of this section.This may be verified by linearizing the adjoint weak form around both the true(which we know for this example) and the approximate solutions. These resultsare shown in Table 5 and now the error estimate is again accurate.11e remark that for the lowest order tuple ( P , P , P ) for the variables( u , b , p ) in Table 2, the error is largely dominated by the contributions E con and E M . We greatly reduce the error in E M by using a higher order space for b as demonstrated in table Table 3. However, this does not reduce the magnitudeof the total error much (about 5%) which is still dominated by the contribution E con . The contribution E con is not significantly affected by the finite dimen-sional space for b . Now finally, if we refine all three variables, Table 4 showsthat the total error drops by two orders of magnitude.2D Elem. True Error Eff. E mom E con E M Table 2: Error in ( u x , Ω c ) for the Hartmann problem § Ω c = [ − , ] × [ − , ].The finite dimensional space here is ( P , P , P ) for ( u , b , p ).
2D Elem. True Error Eff. E mom E con E M Table 3: Error in ( u x , Ω c ) for the Hartmann problem § P , P , P ) for u , b , p ).
2d Elem. True Error Eff. E mom E con E M Table 4: Error in ( u x , Ω c ) for the Hartmann problem § P , P , P ) for ( u , b , p ). Here, we approximate the true solution with the computedsolution which results in linearization error. For this accurate a solution, this deteriorates thequality of the estimate which in turn results in a efficiency further from 1. This is confirmedin Table 5 where we use the true solution and the effectivity is again close to 1. The magnetic lid driven cavity is another common benchmark problem forverfying MHD codes [5, 54]. However, the standard lid velocity is discontinuousand therefore has at most H / − ε regularity in 2D with ε >
0. By the converse12d Elem. True Error Eff. E mom E con E M Table 5: Error in ( u x , Ω c ) for the Hartmann problem, § P , P , P ) for ( u , b , p ). No linearization error is present here because we use the truesoltuion in the definition of the adjoint. of the trace theorem and the Sobelev inequality [55, 56], the solution u x cannotobtain H regularity on the interior. Indeed, in this situation, we do not evenhave well posedness of the primal problem, so there is not real hope for erroranalysis. This issue has been address in a purely fluid context [57, 58]. In bothcases, a regularization of the lid velocity is proposed to mitigate theoreticalissues (in the former) and the ability to achieve higher Reynold’s numbers (inthe latter). In this work, we use a similar regularization to the one proposed in[58], a polynomial regularization of the lid velocity, u top ( x ) = C (cid:0) x − (cid:1) (cid:0) x + (cid:1) , with C chosen such that (cid:90) / − / u top ( x ) d x = 1 . The boundary conditions are imposed as g ( x, .
5) = (cid:2) u top , (cid:3) T on the top faceand zero on the rest of the boundary. The boundary conditions for the magneticfield are q = (cid:2) , (cid:3) T on the left and right ( x = − . x = 0 .
5) and zero onthe top and bottom ( y = − . y = 0 . R f = 5000 and varying magnetic Reynold’s numbers R m in Figure 1. These plots are qualitatively similar to Figure 1 in [5] (forwhich an un-regularized lid velocity is used), which gives us a good indicationthat the regularized version produces qualitatively similar features.Furthermore, since high R f flows provide difficulties for the continuousGalerkin method [59], we use a homotopic sequence of initial guesses to achievehigh R f . Indeed, we run the problem for a moderate value of R f = 200 forexample, and then use the solution produced by the solver as the initial guessfor a larger value e.g. R f = 1000 until we have achieved the desired value.Figure 2 shows the intermediate values in this sequence to solve a problem with R f = 1000. We take the same QoI as for the Hartmann problem although now withΩ c := (cid:2) − , (cid:3) × (cid:2) , (cid:3) , (5.5)13 a) R m = 0 . R m = 0 . R m = 5 . R m . The other nondimensionalized parameters R f = 5000 , κ =1 for all of these plots.(a) R f = 200 (b) R f = 500 (c) R f = 1000Figure 2: Demonstrating the homotopy parameter strategy to achieve high fluid Reynold’snumbers as described in subsection 5.2. The other nondimensionalized parameters R m =5 . , κ = 1 for all of these plots. so that the QoI (Ψ , U ) = ( Ω c , u x ) has the physical interpretation of capturingthe recirculation in the lower half of the box.The effectivity ratio and error contributions various values of R f are shownin Tables 7, 6, 8 and 9. We note that the error estimate is accurate in all caseswith effectivity ratios close to 1.In terms of error contributions, for the lowest order cases, both R f = 2000and R f = 1000 have a fairly balanced distribution of error between the compo-nents as seen in Tables 6 and 8. Similarly to the Hartmann problem, when werefine the space for b , the contribution E M decreases by several orders of mag-nitude, but the overall error remains dominated by E con and E mom . In contrastto Hartmann, E mom is substantial even when we refine b which we attribute tothe high R f nature of the problem. 14d Elem. True Error Eff. E mom E con E M Table 6: error in ( u x , Ω c ) for the lid driven cavity § P , P , P ) for ( u , b , p ). We use an overkill solution on a 400x400=160000 element mesh and( P , P , P ) elements. The parameters are R f = 1000 , R m = 0 . , κ = 1.
2d Elem. True Error Eff. E mom E con E M Table 7: error in ( u x , Ω c ) for the lid driven cavity § P , P ) , P ) for ( u , b , p ). We use an overkill solution on a 400x400=160000 element meshand ( P , P , P ) elements. The parameters are R f = 1000 , R m = 0 . , κ = 1.
2d Elem. True Error Eff. E mom E con E M Table 8: error in ( u x , Ω c ) for the lid driven cavity § P , P ) , P ) for ( u , b , p ). We use an overkill solution on a 400x400=160000 element meshand ( P , P , P ) elements. The parameters are R f = 2000 , R m = 0 . , κ = 1.
2d Elem. True Error Eff. E mom E con E M Table 9: error in ( u x , Ω c ) for the lid driven cavity § P , P , P ) for ( u , b , p ). We use an overkill solution on a 400x400=160000 element meshand ( P , P , P ) elements. The parameters are R f = 2000 , R m = 0 . , κ = 1.
6. Derivation of the weak adjoint and well-posedness
In this section we provide the details of computing the adjoint to exactpenalty weak form following the theory in §
3. Then we use a standard saddlepoint argument to demonstrate the well-posedness of this new adjoint problem(4.5). We take inspiration for these proofs from [3]. To simplify notation in thissection, we define s := u + u h , t := b + b h . (6.1)15inally, we use the notation ( · ) = and ( · ) ≤ to denote that the equality or inequalityis justified by equation ( · ). In this section we provide derivation for the primal linearized operators J ∗ = Y ∗ , J ∗ = Z ∗ u , J ∗ = Z ∗ b and J ∗ = C ∗ in (4.4). We first com-pute the primal linearized operators, Y = J , Z u = J , Z b = J and C = J , using (3.4) and then apply (3.5) to compute the J ∗ ij s. We have from(3.4), Y d := (cid:90) ∂ Y ∂ b ( s b + (1 − s ) b h ) d d s, Z b d := (cid:90) ∂ Z ∂ b ( s u + (1 − s ) u h ) d d s, Z u w := (cid:90) ∂ Z ∂ u ( s b + (1 − s ) b h ) w d s. To this end, we compute Y d = (cid:90) ∂ Y ∂ b ( s b + (1 − s ) b h ) d d s = (cid:90) [ ∇ × ( s b + (1 − s ) b h )] × d + ( ∇ × d ) × ( s b + (1 − s ) b h ) d s = 12 [( ∇ × ( b h + b )) × d + ( ∇ × d ) × ( b h + b )] . (6.2)Similarly, for the two Z terms, Z b d = (cid:90) ∂ Z ∂ b ( s u + (1 − s ) u h ) d d s = (cid:90) ∇ × (( s u + (1 − s ) u h ) × d ) d s = 12 [ ∇ × (( u h + u ) × d )] . (6.3)An identical procedure produces, Z u w = 12 [ ∇ × ( w × ( b + b h ))] . (6.4)Now, to find the adjoints of these operators, we use (3.5), which in our caseinvolves multiplying by a test function and then isolating the trial function usingintegration by parts. We also make use of the vector identities in Appendix B.We are now prepared to compute the adjoint for Y . Integrating (6.2) against16 ∈ H (Ω),( Y d , v ) = 12 (cid:90) Ω [( ∇ × t ) × d + ( ∇ × d ) × t ] · v d x (B.1a) = 12 (cid:90) Ω d · [ v × ( ∇ × t )] + ( ∇ × d ) · [ t × v ] d x (B.1d) = 12 (cid:90) Ω − d · [( ∇ × t ) × v ] + d · [ ∇ × ( t × v )] d x − (cid:90) ∂ Ω d · [( t × v ) × n ] d s (B.1a) = 12 (cid:90) Ω − d · [( ∇ × t ) × v ] + d · [ ∇ × ( t × v )] d x + 12 (cid:90) ∂ Ω ( t × v ) · [ d × n ] d s (2.4) = 12 (cid:90) Ω − d · [( ∇ × t ) × v ] + d · [ ∇ × ( t × v )] d x (4.4) = ( d , Y ∗ v ) . We proceed with computing the adjoint for Z u ,( Z u w , c ) = 12 ( ∇ × ( w × t ) , c ) (B.1d) = 12 (cid:90) Ω ( w × t ) · ( ∇ × c ) d x − (cid:90) ∂ Ω ( w × t ) · ( c × n ) d s (B.1a) = 12 (cid:90) Ω w · [ t × ( ∇ × c )] d x − (cid:90) ∂ Ω ( w × t ) · ( c × n ) d s (2.3) = 12 (cid:90) Ω w · [ t × ( ∇ × c )] d x (4.4) = ( w , Z ∗ u c ) . Finally we compute the adjoint to the linearized operator Z b ,( Z b d , c ) = 12 ( ∇ × ( s × d ) , c ) (B.1d) = 12 (cid:90) Ω ( s × d ) · ( ∇ × c ) d x − (cid:90) ∂ Ω ( s × d ) · ( c × n ) d s (B.1a) = 12 (cid:90) Ω d · [( ∇ × c ) × s ] d x − (cid:90) ∂ Ω d · [ s × ( c × n )] − ( s × d ) · ( c × n ) d s (2.4) = 12 (cid:90) Ω d · [( ∇ × c ) × s ] d x (4.4) = ( d , Z ∗ b c ) . The operator C ∗ is identical to the one presented in [30]. In this section we prove the well-posedness of the adjoint problem § X := H (Ω) × H τ (Ω) and M := L (Ω) so that M (Ω) = X × M . We equip thespace X with the graph norm (cid:107) ( v , c ) (cid:107) X := ( (cid:107) v (cid:107) + (cid:107) c (cid:107) ) / . (6.5)17e next define the bilinear form a : X → R by a (( φ , β ) , ( v , c )) = 1 R f ( ∇ φ , ∇ v ) + (cid:16) C ∗ φ , v (cid:17) + κR m ( ∇ × β , ∇ × c ) + κR m ( ∇ · β , ∇ · c ) − κ (cid:16) Y ∗ φ , c (cid:17) − κ (cid:16) Z ∗ u β , v (cid:17) − κ (cid:16) Z ∗ b β , c (cid:17) , (6.6)and the mixed form b : X × M → R by b (( v , c ) , π ) = − ( π, ∇ · v ) . (6.7)The weak dual problem (4.5) is then equivalent to the following mixed problem:find (( φ , β ) , π ) ∈ X × M such that (cid:40) a (( φ , β ) , ( v , c )) + b (( v , c ) , p ) = ( ψ , v ) , ∀ ( v , c ) ∈ X,b (( φ , β ) , q ) = ( ψ p , q ) , ∀ q ∈ M. (6.8)According the theory of saddle point systems, in order to show the existenceand uniqueness of solutions to (6.8), we must show three things:(i) The bilinear forms a ( · , · ) and b ( · , · ) are bounded on their respective do-mains.(ii) The form a ( · , · ) is coercive on X := { v ∈ X : b ( v, q ) = 0 , ∀ q ∈ M } .(iii) The form b ( · , · ) satisfies the inf-sup condition: ∃ β > q ∈ M sup ( v , c ) ∈ X b (( v , c ) , q ) (cid:107) ( v , c ) (cid:107) X (cid:107) q (cid:107) M ≥ β. (6.9)We organize these parts in the following lemmas. We make frequent use of theinequalities in Appendix C in the proofs. Lemma 1.
The form a ( · , · ) is bounded on X .Proof. Consider the splitting a (( φ , β ) , ( v , c )) = a (( φ , β ) , ( v , c )) + a (( φ , β ) , ( v , c )) (6.10)where a (( φ , β ) , ( v , c )) = 1 R f ( ∇ φ , ∇ v ) + κR m ( ∇ × β , ∇ × c ) + κR m ( ∇ · β , ∇ · c ) ,a (( φ , β ) , ( v , c )) = (cid:16) C ∗ φ , v (cid:17) − κ (cid:16) Y ∗ φ , c (cid:17) − κ (cid:16) Z ∗ u β , v (cid:17) − κ (cid:16) Z ∗ b β , c (cid:17) . Then it suffices to show that both a and a are bounded separately. The prooffor the boundedness of a is given in [3]. For a observe that | a (( φ , β ) , ( v , c )) | ≤ (cid:90) Ω (cid:12)(cid:12)(cid:12) C ∗ φ · v (cid:12)(cid:12)(cid:12) d x + κ (cid:90) Ω (cid:12)(cid:12)(cid:12) Y ∗ φ · c (cid:12)(cid:12)(cid:12) d x + κ (cid:90) Ω (cid:12)(cid:12)(cid:12) Z ∗ u β · v (cid:12)(cid:12)(cid:12) d x + κ (cid:90) Ω (cid:12)(cid:12)(cid:12) Z ∗ b β · c (cid:12)(cid:12)(cid:12) d x. (6.11)18ow, for the first term on the right hand side of (6.11), (cid:90) Ω (cid:12)(cid:12)(cid:12) C ∗ φ · v (cid:12)(cid:12)(cid:12) d x = 12 (cid:90) Ω (cid:12)(cid:12) (cid:2) φ · ∇ s T − ( s · ∇ φ ) − ( ∇ · s ) φ (cid:3) · v (cid:12)(cid:12) d x ≤ (cid:2)(cid:0) φ , ( ∇ s ) T v (cid:1) + ( s , ( ∇ φ ) v ) + ( ∇ · s , φ · v ) (cid:3) ≤ (cid:2) (cid:107) φ (cid:107) L (cid:107) ( ∇ s ) T v (cid:107) L + (cid:107) s (cid:107) L (cid:107) ( ∇ φ ) v (cid:107) L + (cid:107)∇ · s (cid:107) L (cid:107) φ · v (cid:107) L (cid:3) (B.2c) ≤ (cid:2) (cid:107) φ (cid:107) L (cid:107) ( ∇ s ) T (cid:107) L (cid:107) v (cid:107) L + (cid:107) s (cid:107) L (cid:107) ( ∇ φ ) (cid:107) L (cid:107) v (cid:107) L + √ (cid:107)∇ s (cid:107) L (cid:107) φ (cid:107) L (cid:107) v (cid:107) L (cid:3) (C.3) ≤ C Ω (cid:16) (cid:107) φ (cid:107) L (cid:107) s (cid:107) (cid:107) v (cid:107) L + (cid:107) s (cid:107) L (cid:107) φ (cid:107) (cid:107) v (cid:107) L + √ (cid:107) s (cid:107) (cid:107) φ (cid:107) L (cid:107) v (cid:107) L (cid:17) (C.1) ≤ C Ω γ (cid:16) (cid:107) φ (cid:107) (cid:107) s (cid:107) (cid:107) v (cid:107) + (cid:107) s (cid:107) (cid:107) φ (cid:107) (cid:107) v (cid:107) + √ (cid:107) s (cid:107) (cid:107) φ (cid:107) (cid:107) v (cid:107) (cid:17) ≤ √ γC Ω (cid:107) s (cid:107) (cid:107) φ (cid:107) (cid:107) v (cid:107) , where γ is the square of the embedding constant of H (Ω) into L (Ω), see (C.1).For the second term on the right hand side of (6.11), κ (cid:16) Y ∗ φ · c (cid:17) ≤ κ (cid:90) Ω (cid:12)(cid:12) c · [( ∇ × t ) × φ ] (cid:12)(cid:12) + (cid:12)(cid:12) c · [ ∇ × ( t × φ )] (cid:12)(cid:12) d x (B.1d) = κ (cid:90) Ω (cid:12)(cid:12) c · [( ∇ × t ) × φ ] (cid:12)(cid:12) + (cid:12)(cid:12) ( ∇ × c ) · [ t × φ ] (cid:12)(cid:12) d x ≤ κ (cid:107) c (cid:107) L (cid:107) ( ∇ × t ) × φ (cid:107) L + (cid:107)∇ × c (cid:107) L (cid:107) t × φ (cid:107) L ) (C.3) ≤ κC Ω (cid:107) c (cid:107) L (cid:107)∇ × t (cid:107) L (cid:107) φ (cid:107) L + (cid:107)∇ × c (cid:107) L (cid:107) t (cid:107) L (cid:107) φ (cid:107) L ) (B.2d) ≤ κ √ C Ω (cid:107) c (cid:107) L (cid:107) t (cid:107) (cid:107) φ (cid:107) L + (cid:107) c (cid:107) (cid:107) t (cid:107) L (cid:107) φ (cid:107) L ) (C.1) ≤ κγ √ C Ω (cid:107) c (cid:107) (cid:107) t (cid:107) (cid:107) φ (cid:107) . For the third term on the right hand side of (6.11), κ (cid:16) Z ∗ u β , v (cid:17) ≤ κ (cid:90) Ω (cid:12)(cid:12) v · [ t × ( ∇ × β )] (cid:12)(cid:12) d x ≤ κ (cid:107) v (cid:107) L (cid:107) t × ( ∇ × β ) (cid:107) L ≤ κ (cid:107) v (cid:107) L (cid:107) t (cid:107) L (cid:107)∇ × β (cid:107) L (C.3) , (B.2d) ≤ κC Ω √ (cid:107) v (cid:107) L (cid:107) t (cid:107) L (cid:107) β (cid:107) ≤ κC Ω γ √ (cid:107) v (cid:107) (cid:107) t (cid:107) (cid:107) β (cid:107) . The fourth term follows the same argument as the third term to yield the bound, κ (cid:16) Z ∗ b β , c (cid:17) ≤ κC Ω γ √ (cid:107) c (cid:107) (cid:107) s (cid:107) (cid:107) β (cid:107) . (6.12)19utting these bounds together, we conclude a (( φ , β ) , ( v , c )) ≤ γC Ω (cid:18) √ (cid:107) s (cid:107) (cid:107) φ (cid:107) (cid:107) v (cid:107) + κ √ (cid:107) c (cid:107) (cid:107) t (cid:107) (cid:107) φ (cid:107) + κ √ (cid:107) v (cid:107) (cid:107) t (cid:107) (cid:107) β (cid:107) + κ √ (cid:107) c (cid:107) (cid:107) s (cid:107) (cid:107) β (cid:107) (cid:19) (C.5) ≤ γC Ω (cid:32) √ (cid:107) s (cid:107) (cid:107) φ (cid:107) (cid:107) v (cid:107) + κ √ (cid:107) c (cid:107) (cid:107) s (cid:107) (cid:107) β (cid:107) + (cid:107) t (cid:107) κ √ (cid:107) ( v , c ) (cid:107) X (cid:107) ( φ , β ) (cid:107) X (cid:33) (C.5) ≤ γC Ω (cid:32) (cid:107) s (cid:107) max (cid:40) √ , κ √ (cid:41) (cid:107) ( v , c ) (cid:107) X (cid:107) ( φ , β ) (cid:107) X + (cid:107) t (cid:107) (cid:107) ( v , c ) (cid:107) X (cid:107) ( φ , β ) (cid:107) X (cid:33) ≤ α b (cid:107) ( v , c ) (cid:107) X (cid:107) ( φ , β ) (cid:107) X , (6.13)where in turn α b = max (cid:40) (cid:107) s (cid:107) max (cid:40) √ , κ √ (cid:41) , (cid:107) t (cid:107) (cid:41) . Now we consider the coercivity of the bilinear form a ( · , · ) on X . Lemma 2.
There exists a constant α c > such that whenever k R f − γC Ω (cid:34) √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35) > , (6.14) and k κR m − γC Ω (cid:34) κ √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35) > then a (( φ , β ) , ( φ , β )) ≥ α c (cid:107) ( φ , β ) (cid:107) X , ∀ ( φ , β ) ∈ X. (6.16) Proof.
Using the splitting established in the previous lemma, a (( φ , β ) , ( φ , β )) ≥ a (( φ , β ) , ( φ , β )) − | a (( φ , β ) , ( φ , β )) | = 1 R f ( ∇ φ , ∇ φ ) + κR m ( ∇ × β , ∇ × β ) + κR m ( ∇ · β , ∇ · β ) − | a (( φ , β ) , ( φ , β )) |≥ k R f (cid:107) φ (cid:107) + k κR m (cid:107) β (cid:107) − | a (( φ , β ) , ( φ , β )) | (6.17)20here k comes from the Poincar´e inequality(C.6), and k is defined though (cid:107)∇ × v (cid:107) + (cid:107)∇ · v (cid:107) ≥ k (cid:107) v (cid:107) , ∀ v ∈ H τ (Ω) , (6.18)which is valid under the restrictions we have imposed on the domain Ω and thecontinuous embed ding of H τ (Ω) (cid:44) → H (Ω) [42]. Picking up from (6.17) andusing (C.7) we conclude that, a (( φ , β ) , ( φ , β )) ≥ k R f (cid:107) φ (cid:107) + k κR m (cid:107) β (cid:107) ≥ (cid:32) k R f − γC Ω √ (cid:107) s (cid:107) (cid:33) (cid:107) φ (cid:107) + (cid:32) k κR m − γC Ω κ √ (cid:107) s (cid:107) (cid:33) (cid:107) β (cid:107) − γC Ω κ √ (cid:107) φ (cid:107) (cid:107) t (cid:107) (cid:107) β (cid:107) ≥ (cid:32) k R f − γC Ω √ (cid:107) s (cid:107) (cid:33) (cid:107) φ (cid:107) + (cid:32) k κR m − γC Ω κ √ (cid:107) s (cid:107) (cid:33) (cid:107) β (cid:107) − γC Ω κ √ (cid:107) t (cid:107) (cid:0) (cid:107) β (cid:107) + (cid:107) φ (cid:107) (cid:1) = (cid:32) k R f − γC Ω (cid:34) √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35)(cid:33) (cid:107) φ (cid:107) + (cid:32) k κR m − γC Ω (cid:34) κ √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35)(cid:33) (cid:107) β (cid:107) . Thus, taking α c = min (cid:40) k R f − γC Ω (cid:34) √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35) ,k κR m − γC Ω (cid:34) κ √ (cid:107) s (cid:107) + 3 κ √ (cid:107) t (cid:107) (cid:35) (cid:41) , (6.19)concludes the lemma.Now we are prepared to prove the main result. Theorem 3.
Under the conditions of Lemma 2 there exists a unique solutionto the dual problem (4.5) .Proof.
The boundedness and inf-sup condition for b ( · , · ) are standard see e.g.[56]. The boundedness of a ( · , · ) follows from Lemma 1, and Lemma 2 proves a ( · , · ) is coercive on X so in particular on X .21 ppendix A. Standard function spaces We denote by L (Ω) the set of all square Lebesgue integrable functions onΩ ⊂ R d with associated inner product ( · , · ) and norm (cid:107)·(cid:107) . This extends naturallyto vector valued functions, denoted by L (Ω), where the inner product is givenby, ( u , v ) = d (cid:88) i =1 ( u i , v i ) . (A.1)The Sobolev norm for p = 2 is, (cid:107) v (cid:107) m := m (cid:88) | α | =0 (cid:13)(cid:13) D α v (cid:13)(cid:13) / . where α = ( α , . . . , α m ) is a multi-index of length m and D α v := ∂ α x . . . ∂ α m x m v, where the partial derivatives are taken in the weak sense. Thus, the Hilbertspaces H m for m = 0 , , , . . . is simply be defined as functions with bounded m -norm, H m (Ω) := { v : (cid:107) v (cid:107) m < ∞} . The space H (Ω) is identified with L (Ω). For vector valued functions, theHilbert space H m is defined as, H m (Ω) := { v : v i ∈ H m (Ω) , i = 1 , . . . , d } . Appendix B. Vector identities and inequalities
In this section we have collected all relevant identities to perform the neces-sary integration by parts arguments for the readers convenience, A · ( B × C ) = B · ( C × A ) = C · ( A × B ) , (B.1a) ∇ · ( A × B ) = B · ( ∇ × A ) − A · ( ∇ × B ) , (B.1b) ∇ · [( ∇ · A ) B ] = ( ∇ · A )( ∇ · B ) + B · [ ∇ ( ∇ · A )] , (B.1c) (cid:90) Ω A · ( ∇ × B ) d x = (cid:90) Ω B · ( ∇ × A ) d x − (cid:90) ∂ Ω B · ( A × n ) d s, (B.1d) (cid:90) Ω B · [ ∇ ( ∇ · A )] d x = − (cid:90) Ω ( ∇ · A )( ∇ · B ) d x + (cid:90) ∂ Ω ( ∇ · A ) B d s. (B.1e)One should note that the integral identities (B.1d) and (B.1e) follow from thecomponent-wise identities (B.1a)-(B.1c) and the divergence theorem. We also22ake use of the following inequalities for u , v ∈ H , | u · v | ≤ (cid:107) u (cid:107) R (cid:107) v (cid:107) R , (B.2a) (cid:107) u × v (cid:107) ≤ (cid:107) u (cid:107) R (cid:107) v (cid:107) R , (B.2b) (cid:107)∇ × u (cid:107) ≤ √ (cid:107)∇ u (cid:107) R × R , (B.2c) (cid:107)∇ · u (cid:107) ≤ √ (cid:107)∇ u (cid:107) R × R (B.2d) (cid:107) Av (cid:107) R ≤ (cid:107) A (cid:107) R × (cid:107) v (cid:107) R , (B.2e)and finally the equality (cid:107)∇ v T (cid:107) R × = (cid:107)∇ v (cid:107) R × . (B.3) Appendix C. Useful inequalities from analysis
1. The space H (Ω) embeds continuously in L (Ω) with constant √ γ . Thatis, H (Ω) (cid:44) → L (Ω) such that, (cid:107) v (cid:107) L ≤ √ γ (cid:107) v (cid:107) H . (C.1)2. We have (cid:107) u · v (cid:107) L ≤ (cid:107) u (cid:107) L (cid:107) v (cid:107) L . (C.2)This is seen as follow, (cid:107) u · v (cid:107) L = (cid:18)(cid:90) Ω ( u · v ) d x (cid:19) / ≤ (cid:18)(cid:90) Ω (cid:107) u (cid:107) R (cid:107) v (cid:107) R d x (cid:19) / ≤ (cid:32)(cid:18)(cid:90) Ω (cid:107) u (cid:107) R × (cid:19) / (cid:18)(cid:90) Ω (cid:107) v (cid:107) R (cid:19) / (cid:33) / = (cid:107) u (cid:107) L (cid:107) v (cid:107) L .
3. For v ∈ L , we have (cid:107) v (cid:107) L ≤ C Ω (cid:107) v (cid:107) L . (C.3)Here C Ω is a constant that depends on Ω. The last inequality follows fromthe H¨older inequality. Suppose r satisfies q + r = p . From the H¨olderinequality, (cid:107) f (cid:107) pp = (cid:107) f (cid:107) pp = (cid:107) p f p (cid:107) ≤ (cid:107) f p (cid:107) q/p (cid:107) p (cid:107) r/p = (cid:107) f (cid:107) pq (cid:107) (cid:107) pr = ⇒ (cid:107) f (cid:107) p ≤ meas(Ω) r (cid:107) f (cid:107) q . In this specific case, q = r = 4 and p = 2.4. We have, (cid:107) u × v (cid:107) L ≤ √ (cid:107) u (cid:107) L (cid:107) v (cid:107) L (C.4)23s shown below, (cid:107) u × v (cid:107) L = (cid:18)(cid:90) Ω (cid:107) u × v (cid:107) R d x (cid:19) / ≤ (cid:18)(cid:90) Ω √ (cid:107) u (cid:107) R (cid:107) v (cid:107) R d x (cid:19) / ≤ (cid:32) √ (cid:18)(cid:90) Ω (cid:107) u (cid:107) R × (cid:19) / (cid:18)(cid:90) Ω (cid:107) v (cid:107) R (cid:19) / (cid:33) / = √ (cid:107) u (cid:107) L (cid:107) v (cid:107) L .
5. The Cauchy-Schwarz inequality for (cid:2) a, b (cid:3) , (cid:2) c, d (cid:3) ∈ R , ac + bd = (cid:2) a, b (cid:3) (cid:2) c, d (cid:3) T ≤ (cid:112) a + c (cid:112) b + d , (C.5)6. The Poincar´e inequality is, (cid:107)∇ v (cid:107) ≥ k (cid:107) v (cid:107) , ∀ v ∈ H (Ω) . (C.6)7. For x, y ∈ R , − xy ≥ − ( x + y ) , (C.7) Acknowledgements
J. Chaudhry’s work is supported by the NSF-DMS 1720402. The work ofAri. E. Rappaport and John N. Shadid was partially supported by the U.S.Department of Energy, Office of Science, Office of Advanced Scientific Com-puting Research, Applied Mathematics Program and by the U.S. Departmentof Energy, Office of Science, Office of Advanced Scientific Computing Researchand Office of Fusion Energy Sciences, Scientific Discovery through AdvancedComputing (SciDAC) program. Sandia National Laboratories is a multimissionlaboratory managed and operated by National Technology and Engineering So-lutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International,Inc., for the U.S. Department of Energys National Nuclear Security Administra-tion under contract DE-NA0003525. This paper describes objective technicalresults and analysis. Any subjective views or opinions that might be expressedin the paper do not necessarily represent the views of the U.S. Department ofEnergy or the United States Government.
References [1] J. P. H. Goedbloed, S. Poedts, A. Mills, S. Romaine, Principles of Mag-netohydrodynamics: With Applications to Laboratory and AstrophysicalPlasmas, Cambridge University Press, 2003.[2] M. Ulrich, B. Leo, Magnetofluiddynamics in Channels and containers,Springer, 2010. 243] M. D. Gunzburger, A. J. Meir, J. S. Peterson, On the existence, uniqueness,and finite element approximation of solutions of the equations of station-ary, incompressible magnetohydrodynamics, Mathematics of Computation56 (194) (1991) 523–523.[4] J. Shadid, R. Pawlowski, J. Banks, L. Chac´on, P. Lin, R. Tuminaro, To-wards a scalable fully-implicit fully-coupled resistive MHD formulation withstabilized FE methods, Journal of Computational Physics 229 (20) (2010)76497671.[5] E. G. Phillips, H. C. Elman, E. C. Cyr, J. N. Shadid, R. P. Pawlowski, Ablock preconditioner for an exact penalty formulation for stationary MHD,SIAM Journal on Scientific Computing 36 (6) (2014).[6] P. T. Lin, J. N. Shadid, P. H. Tsuji, On the performance of krylov smooth-ing for fully coupled AMG preconditioners for VMS resistive MHD, Inter-national Journal for Numerical Methods in Engineering 120 (12) (2019)1297–1309.[7] P. Lin, J. Shadid, J. Hu, R. Pawlowski, E. Cyr, Performance of fully-coupled algebraic multigrid preconditioners for large-scale VMS resistiveMHD, Journal of Computational and Applied Mathematics 344 (2018) 782–793.[8] J. Shadid, R. Pawlowski, J. Banks, L. Chac´on, P. Lin, R. Tuminaro, To-wards a scalable fully-implicit fully-coupled resistive MHD formulation withstabilized fe methods, Journal of Computational Physics 229 (20) (2010)7649–7671.[9] L. Chac´on, A. Stanier, A scalable, fully implicit algorithm for the reducedtwo-field low- β extended mhd model, Journal of Computational Physics326 (2016) 763–772.[10] E. C. Cyr, J. N. Shadid, R. S. Tuminaro, R. P. Pawlowski, L. Chac´on,A new approximate block factorization preconditioner for two-dimensionalincompressible (reduced) resistive MHD, SIAM Journal on Scientific Com-puting 35 (3) (2013) B701–B730.[11] J. H. Adler, T. A. Manteuffel, S. F. McCormick, J. W. Ruge, First-ordersystem least squares for incompressible resistive magnetohydrodynamics,SIAM Journal on Scientific Computing 32 (1) (2010) 229–248.[12] J. H. Adler, T. A. Manteuffel, S. F. McCormick, J. W. Ruge, G. D. Sanders,Nested iteration and first-order system least squares for incompressible,resistive magnetohydrodynamics, SIAM Journal on Scientific Computing32 (3) (2010) 1506–1526.[13] J. H. Adler, M. Brezina, T. A. Manteuffel, S. F. McCormick, J. W. Ruge,L. Tang, Island coalescence using parallel first-order system least squares25n incompressible resistive magnetohydrodynamics, SIAM Journal on Sci-entific Computing 35 (5) (2013) S171–S191.[14] P.-W. Hsieh, S.-Y. Yang, A bubble-stabilized least-squares finite elementmethod for steady MHD duct flow problems at high hartmann numbers,Journal of Computational Physics 228 (22) (2009) 8301–8320.[15] C. Trenchea, Unconditional stability of a partitioned IMEX method formagnetohydrodynamic flows, Applied Mathematics Letters 27 (2014) 97–100.[16] R. Codina, N. H. Silva, Stabilized finite element approximation of thestationary magneto-hydrodynamics equations, Computational Mechanics38 (4-5) (2006) 344–355.[17] J. Li, H. Yang, E. Machorro (Eds.), Recent Advances in Scientific Comput-ing and Applications, American Mathematical Society, 2013.[18] D. Estep, V. Ginting, D. Ropp, J. N. Shadid, S. Tavener, An a posteriori–apriori analysis of multiscale operator splitting, SIAM Journal on NumericalAnalysis 46 (3) (2008) 1116–1146.[19] D. Estep, A. M˚alqvist, S. Tavener, Nonparametric density estimation forrandomly perturbed elliptic problems I: Computational methods, a poste-riori analysis, and adaptive error control, SIAM J. Sci. Comput. 31 (2009)2935–2959.[20] J. H. Chaudhry, N. Burch, D. Estep, Efficient distribution estimation anduncertainty quantification for elliptic problems on domains with stochasticboundaries, SIAM/ASA Journal on Uncertainty Quantification 6 (3) (2018)1127–1150.[21] D. Estep, A posteriori error bounds and global error control for approxima-tion of ordinary differential equations, SIAM Journal on Numerical Analysis32 (1) (1995) 1–48.[22] M. B. Giles, E. Sli, Adjoint methods for PDEs: a posteriori error analysisand postprocessing by duality, Acta Numerica 2002 (2002) 145236.[23] M. Ainsworth, T. Oden, A posteriori error estimation in finite elementanalysis, John Wiley-Teubner, 2000.[24] W. Bangerth, R. Rannacher, Adaptive Finite Element Methods for Differ-ential Equations, Birkhauser Verlag, 2003.[25] D. J. Estep, M. G. Larson, R. D. Williams, A. M. Society, Estimatingthe error of numerical solutions of systems of reaction-diffusion equations,American Mathematical Society, 2000.[26] R. Becker, R. Rannacher, An optimal control approach to a posteriori errorestimation in finite element methods: Acta numerica (Jan 2003).2627] Y. Cao, L. Petzold, A posteriori error estimation and global error controlfor ordinary differential equations by the adjoint method, SIAM Journalon Scientific Computing 26 (2) (2004) 359–374.[28] J. H. Chaudhry, J. N. Shadid, T. Wildey, A posteriori analysis of an IMEXentropy-viscosity formulation for hyperbolic conservation laws with dissi-pation, Applied Numerical Mathematics 135 (2019) 129–142.[29] J. H. Chaudhry, J. Collins, J. N. Shadid, A posteriori error estimation formulti-stage runge–kutta IMEX schemes, Applied Numerical Mathematics117 (2017) 36–49.[30] D. Estep, S. Tavener, T. Wildey, A posteriori error estimation and adaptivemesh refinement for a multiscale operator decomposition approach to fluid-solid heat transfer, Journal of Computational Physics 229 (11) (2010) 4143–4158.[31] K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Introduction to adaptivemethods for differential equations, Acta Numerica 4 (1995) 105–158.[32] J. H. Chaudhry, A posteriori analysis and efficient refinement strategies forthe poisson–boltzmann equation, SIAM Journal on Scientific Computing40 (4) (2018) A2519–A2542.[33] K. Eriksson, D. Estep, P. Hansbo, C. Johnson, Computational DifferentialEquations, Cambridge University Press, Cambridge, 1996.[34] T. J. Barth,
A posteriori
Error Estimation and Mesh Adaptivity for FiniteVolume and Finite Element Methods, Vol. 41 of Lecture Notes in Compu-tational Science and Engineering, Springer, New York, 2004.[35] J. H. Chaudhry, D. Estep, V. Ginting, J. N. Shadid, S. Tavener, A posteriorierror analysis of imex multi-step time integration methods for advection–diffusion–reaction equations, Computer Methods in Applied Mechanics andEngineering 285 (2015) 730–751.[36] D. Schotzau, Mixed finite element methods for stationary incompressiblemagneto-hydrodynamics, Numerische Mathematik 96 (4) (2004) 771800.[37] J. C. Nedelec, Mixed finite elements in R3