Constraint preserving boundary conditions for the Ideal Newtonian MHD equations
aa r X i v : . [ a s t r o - ph ] J a n Constraint preserving boundary conditions for the Ideal NewtonianMHD equations
Mariana C´ecere , Luis Lehner and Oscar Reula FaMAF, Universidad Nacional de C´ordoba, C´ordoba, 5000 (Argentina) Department of Physics and Astronomy, Louisiana State University, Baton Rouge, LA 70803-4001
1. Introduction
Magnetic fields play an important role in the behavior of plasmas and are thought mediateimportant effects like dynamos in the core of planets and the formation of jets in active galacticnuclei and gamma ray bursts; induce a variety of magnetic instabilities; realize solar flares, etc.(see e.g. [1,2]). Understanding the role of magnetic fields in these and other phenomena havespurred through the years many efforts to obtain solutions of the magnetohydrodynamic equa-tions. The non-linear nature of these equations limits the understanding that can be gainedin a particular problem via analytical techniques. This implies that solutions for complex sys-tems must be obtained by numerical means and a suitable numerical implementation must beconstructed for this purpose. Such implementation must be able to evolve the solution to thefuture of some initial configuration and guarantee its quality. A delicate, subsidiary quantity,can be monitored in part to estimate this. This quantity is the “monopole constraint” ∂ i B i which must be zero at the analytical level for a consistent solution. This quantity is not a partof the main variables, rather it is a derived quantity which should be satisfied by a true physi-cal solution. In practice, unless a numerical implementation of the MHD equations is carefullydesigned this constraint can be severely violated. This, in turn, signals (and is sometimes thecause) of a degrading numerical solution. For these reasons, several approaches have been in-vestigated and developed for guaranteeing a controlled behavior of this quantity. One suchapproach is known as the constraint transport technique [3,4] which adopts a particular algo-rithm that staggers the variables appropriately to ensure the satisfaction of the constraint atround-off level within Finite Difference and Finite Elements techniques. This approach hasbeen quite successful in a number of applications across different disciplines and particularlyrelevant in astrophysics applications [5,6,7,8,9,10,11,12]. However, by design it imposes limitson the algorithmic options available to an implementation. This fact can be at odds, or intro-duce complications, with applications where adaptive mesh refinement is required and/or ad-vanced numerical techniques (that exploit useful properties of the equations) are adopted. Analternative approach, which controls the constraint at truncation-error level maintains com-plete freedom in the numerical techniques to be adopted. This approach, referred to as diver-gence cleaning puts the burden to control the constraint not on the algorithm to be employed ut rather on the system of equations to be solved itself [13,14]. This is achieved by consider-ing an additional variable suitably coupled to the system through another equation so that,through the evolution, the constraint behavior is kept under control. While this method has,to date, received less attention than the constraint transport one, successful applications indiverse scenarios have already illustrated its usefulness (e.g. [14,15,4,16,17]).Regardless of the technique employed, boundary conditions can play a significant role andmay spoil all efforts towards a correct implementation in the absence of boundaries. Clearly,even when a stable option is adopted, it need not be consistent with the goal of preservingthe constraint (be it at round-off or truncation levels) and so carefully designed boundaryconditions must be formulated. Commonly employed options include straightforward outflow-type conditions –which do not enforce the constraint– or “absorbing” boundary conditionswhich aim to reduce the influence of spurious effects induced at the boundary in the numericalsolution [18].In the present work we concentrate on formulating constraint preserving boundary condi-tions for the Newtonian ideal MHD equations (for systems with and without divergence clean-ing). Such a task is intimately related to the hyperbolic properties of the equations and sowe re-analyze the system of equations and discuss alternatives for controlling the constraintsthrough the divergence cleaning technique.We therefore organize our presentation along the following lines. In section II we analyze thesystem of equations and the formulation of boundary conditions beginning with a simplifiedmodel from which we draw the strategy to apply in the complete MHD system. Section IIIpresents a series of tests that highlight the benefits gained by our construction. We concludein section IV with some final comments.
2. The Newtonian equations & constraint hyperbolic cleaning
The equations describing the ideal Newtonian MHD equations in terms of the variables( ρ, p, v i , B i ) are [1,2]): ∂ t ρ = −∇ i ( ρv i ) ,∂ t p = − v i ∇ i p − γp ∇ j v j ,ρ∂ t v i = − ρv j ∇ j v i − ∇ i p − ∇ j ( e ij B k B k − B j B i ) ,∂ t B i = −∇ j ( v j B i − v i B j ) . In this form, this system of equations is weakly hyperbolic since there is no complete set ofeigenvectors. This indicates instabilities are likely to arise unless some modes are carefullycontrolled. The constraint transport technique attempts to do so at the algorithmical level byenforcing ∂ i B i = 0 at the discrete level. We here choose an alternative approach where thesituation is remedied at the analytical level through the addition of an extra field coupled tothe system in a suitable way [13,14]. To simplify the discussion, we first consider a simplermodel which captures essential features of the main system. Preliminaries
To illustrate this approach consider first the system obtained when the fluid field variablesas given and stationary. The resulting system describes a simple model that shares key prob-lematic features of the original ideal MHD system (1), although in the latter case additionalaspects enter under consideration. In this simple case one just has the magnetic field whichevolves under, ∂ t B i = −∇ j ( v j B i − v i B j ) . (1)This equation is already weakly hyperbolic, since a plane-wave analysis indicates that whenthe wave number vector is perpendicular to the velocity vector there exists only two linearlyindependent eigenvectors. To see this consider the case where the wave number vector is normalto the velocity vector, and define coordinates axis such that the first axis lies along the wavevector –so k i = ( k, , v i = (0 , v, B i = U e σt + k i x i , gives rise to the following problem: σU = M U (with σ and U afrequency and a vector to be determined) and M defined as M = kv (2)which is clearly non diagonalizable. Thus, there is no stability in the usual L norm sense.Interestingly the divergence of the magnetic field is preserved by equation (1), that is, if onedefines D := ∇ k B k and takes the divergence of (1) one obtains the induced evolution equation for D as ∂ t D = 0. Hence, if D is zero initially it remains so as along as the time integrationlines do not intersect a boundary. However, we note that if instead of considering the L norm,one considers a different norm defined by adding a term proportional to the divergence of B i ,namely a norm of the type: E := Z Σ t [ B + c ( ∇ k B k ) ] dV one can show that ˙ E = 0 and so this energy is controlled. Therefore, from an analytical pointof view equation (1) is not a bad one.At the numerical level, things become more delicate as generic violations of this constraintwill arise due to truncation or round-off errors which might grow unless a careful implemen-tation of the equations is adopted. Furthermore, even when an integration scheme is availablethat controls the constraint in the absence of boundaries, (which is not difficult to obtain in flatspace in Cartesian coordinates –which we use in this work–), the boundary conditions modifythe equations at the boundaries which can cause the constraint to grow and propagate to theinterior. Thus, we next discuss a way to achieve a more robust behavior. Rather than doing soby the particular numerical algorithm employed we work first at the analytical level and adopta strategy that would help even beyond the simple problem adopted here. o remedy the lack of strong hyperbolicity and good constraint propagation a possible ap-proach can be adopted which makes use of the freedom to add to equation (1) any term pro-portional to that divergence. The proposed system is (in analogy with that in [14]): ∂ t B i = −∇ j ( v j B i − v i B j ) − c l ∇ i φ ,∂ t φ = − c l ∇ k B k − sφ . (3)Notice that for φ = 0 the modified system is equivalent to the original one, and φ hasas a source proportional to the constraint. Thus if boundary and initial data are such as toguarantee that ∇ k B k = 0, and trivial data is given for φ the solution to the modified systemis also a solution of the original equation. The advantage of this modification is that now thesystem is strongly hyperbolic, so it has a well posed initial problem and is stable irrespectiveof whether or not ∇ k B k vanishes. Furthermore, the induced evolution equation for D is nolonger trivial, and it implies that the field D propagates with speed c l . Consequently one canmake the constraint violations to propagate away from the integration region, and, if correctboundary conditions are given, leave the computational region entirely. Additionally, the lastterm in the equation for φ in (3) induces a decay of the constraint for s > ∂ t B i = − u j ∇ j B i + B j ∇ j u i − B i ∇ j u j + (1 − α ) u i ∇ j B j − c l ∇ i φ∂ t φ = − βu j ∇ j φ − c l ∇ i B i − sφ . (4)Which, for a range of values of the parameters α and β for is strongly hyperbolic. Amongthese, of particular interest is the system with α = β = 1, which is strongly hyperbolic evenin the limit c l →
0, so its associated initial value problem would be well posed even withoutthe coupling to the field φ . Another interesting choice is the one α = β = 0 since the resultingsystem can be expressed in conservative form. This property is often preferred in applicationswhere shocks are present as one can easily take advantage of special algorithms defined to dealwith these issues (see e.g. [19,20,21,22,23]).Boundary conditions To fix ideas, we now study the possible boundary conditions for theabove system (with α = β = 0). Our goal is to define conditions which can yield a wellposed problem and, if possible, guarantee no violations are introduced into the computationaldomain. To this end we must determine which are the incoming and outgoing modes off theboundary surface as this information is key to understand what data is freely specifiable. Forthis we seek a solution of the form ( B i , φ ) = U e σt + n i x i , where U is a vector to be determinedalong with the frequency σ , and n i is the outgoing normal to the boundary under consideration.Applying this solution to the system we obtain: σB i = − B i v n + v i B n − c l n i φ ,σφ = − c l B n , (5) here B n := B i n i , and v n := v i n i . This problem of eigenvalues/eigenvectors can be solved todetermine the incoming ( σ > σ = 0) and outgoing ( σ <
0) modes. A necessarycondition for a well posed problem for hyperbolic systems indicates that data must be givenonly to the incoming modes, since the others are determined from the inside of the integrationregion [25,24]. Recall that the incoming modes can be defined even as linear functions of theoutgoing ones so long as the coefficients are small enough.At a given boundary, there are three eigenvalues, ( σ := − v n , σ ± := ± c l ). The first corre-sponds to two linearly independent modes which span the tangent space of the boundary point(i.e. they are perpendicular to n i ). They are positive whenever the velocity is incoming, theseare the modes which are dragged along v i . The other eigenvalues correspond to the cases wherewe can choose B n := b arbitrarily and, φ = ∓ b ,B i = ( v i ± n i c l ) bv n ± c l . Notice one can always adopt a value for c l larger than the maximum velocity expected so thatdenominator is never zero. The eigenvector corresponding to the positive eigenvalue determinethe combination of variable that must be always suitably given to preserve the constraint. Thecomplete expression for the eigenbase and its co-base is. U := (0 , e , , σ = − v n ,U := (0 , e , , σ = − v n ,U ± := (1 , ( v i ± n i c l ) v n ± c l , ∓ , σ = ± c l , where the generic vector is U := ( B n , ˜ B i , φ ) with ˜ B i n i = 0, and { e , e } are two orthonormalvectors on the tangent of the boundary.The co-base is given by θ := ( v n v c l − v n , e , c l v c l − v n ) ,θ := ( v n v c l − v n , e , c l v c l − v n ) ,θ ± := 12 (1 , , ∓ , where v := e i v i , and v := e i v i . In order to see how we must define boundary values consistentwith the constraints we now analyze the subsidiary system determining the evolution of theconstraint. In order to treat it as a first order system we define a new variable δ i := ∇ i φ , whoseevolution equation is determined by (the time derivative of) the evolution equation of φ . Thefull subsidiary system is then: ∂ t D = − c l ∇ i δ i ,∂ t δ i = − c l ∇ i D − sδ i . (6)Here again we must investigate now which are the incoming and outgoing modes for this system.If we can impose boundary conditions to the main system so that we ensure no incoming mode s created for this subsidiary system, then uniqueness of the (trivial) solution would guaranteethat the solution of the main system has vanishing constraint quantities and so it is a solution ofthe original system. As done previously, plane wave solutions of the above system perpendicularto the boundary give rise to the following eigenvalue/eigenvector problem: κD = − c l n i δ i ,κδ i = − c l n i D , (7)and its eigenvalues are { κ = 0 , κ ± = ± c l } with the corresponding eigenvectors: V := (0 , e , , κ = 0 ,V := (0 , e , , κ = 0 ,V ± := (1 , ( v i ± n i c l ) v n ± c l , ∓ , κ ± = ± c l , where the generic vector is V := ( δ n , ˜ δ i , D ) with ˜ δ i n i = 0, and { e , e } are two orthonormalvectors tangent to the boundary.The co-base is given byΘ := (0 , e , , Θ := (0 , e , , Θ ± := 12 (1 , , ∓ . Here again the theory of boundary conditions indicate we should give as a boundary conditionthe incoming mode as a linear function of the outgoing one. That is we should ensure that theboundary data for the main system is such that:Θ + ( V ) + a Θ − ( V ) ˆ=0 , | a | < D and δ n in terms of the time derivatives of B n and φ (and the remaining equations).Hence δ n = − c l ( ∂ t B n − F n ) ,D = − c l ( ∂ t φ − F φ ) , (9)with F n = − n i ∂ j ( u j B i − u i B j ), and F φ = − sφ .We can then translate the above boundary condition into: − ˙ B n (1 + a ) + ˙ φ (1 − a ) + F n (1 + a ) − F φ (1 − a ) ˆ= 0 , (10)which can be re-expressed as: − θ + ( U ) + aθ − ( U )) + F n (1 + a ) − F φ (1 − a ) ˆ= 0 . (11) or the particular case of setting the incoming constraint modes to zero ( a = 0) this implies, θ + ( U ) ˆ= 12 ( F φ − F n ) . (12)Thus, the boundary condition defined by the equation above is not only maximally dissipativebut also enforces the constraint. The complete system
We now turn our attention to the complete system and develop an analogous strategy fol-lowing the discussion considered above. In most cases the literature on the subject considersthe MHD system expressed in terms variables U = ( ρ, e, u i , B i ) rather than U = ( ρ, p, u i , B i )–ie the internal energy e instead of the pressure p –. While this involves a simple change of vari-ables it turns out the characteristic decomposition is simpler in the latter case. So, our evolu-tion system will be given by the former, and we will obtain the characteristic structure for thelatter followed by a straightforward change of variables. Our system of interest is a variationof the initial MHD equations as (see [14]): ∂ t ρ = −∇ i ( ρv i ) ,ρ∂ t v i = − ρv j ∇ j v i − ∇ i p − B k ( ∇ i B k − ∇ k B i ) − αB i ∇ k B k ,∂ t B i = −∇ j ( u j B i − u i B j ) − αu i ∇ j B j − c l ∇ i φ ,∂ t e = −∇ i (( e + p + 12 B ) v i − B i v · B ) − αv i B i ∇ k B k − c l B k ∇ k φ ,∂ t φ = − αu j ∇ j φ − c l ∇ j B j − sφ , (13)where p := ( γ − (cid:18) e − ρv − B (cid:19) , c s := γpρ (14)Thys system includes both the possibility of divergence cleaning ( c l = 0) and Galilean invari-ance ( α = 1) and is strongly hyperbolic, hence has a complete set of eigenvectors. This willallow us to introduce a new boundary treatment which ensures no constraint violations areintroduced through the computational domain boundaries. The parameter α controls the free-dom of adding the constraint equation. In what follows, we shall use the values α = 0 ,
1. Thevalue α = 0 corresponds to the purely conservative system (notice that velocity equation canbe re-written with the help of the first equation as momentum conservation). The case with α = 1 gives rise to a system that is Galilean invariant (see [14]) and also to the case where thesystem is strongly hyperbolic [26,27] irrespective of the field φ . Hence one can study the limitwhere the φ field decouples and the constraint propagates along the velocity field.Boundary system As done previously, we now consider a boundary with outgoing unit normal n i and obtain the characteristic decomposition at this boundary. For notational purposes wewill employ an overbar to denote the perturbation on a given variable, for instance ¯ ρ will enote the perturbation of ρ which will be considered as a fixed background quantity. Thecharacteristic decomposition is then determined from the system, σ ¯ ρ = − v n ¯ ρ − ρ ¯ v n ,σ ¯ v i = − v n ¯ v i + B i ¯ B n /ρ + ( B n /ρ ) ¯ B i − n i ( B j ¯ B j + ¯ p ) /ρ ,σ ¯ B i = − v n ¯ B i + v i ¯ B n − B i ¯ v n + B n ¯ v i + n i ¯ φ ,σ ¯ p = − c s ρ ¯ v n − ( γ − B j v j ¯ B n − v n ¯ p + v (1 − γ ) B n ¯ φ .σ ¯ φ = c l ¯ B n (15)The solution to this system is cumbersome and requires dealing with different cases. The fullanalysis and solution is presented in appendix A and expressed in terms of a vector whose en-tries are: U = ( ¯ ρ, ¯ v n , ˜¯ v i , ¯ B n , ˜¯ B i , ¯ p, ¯ φ ) T indicating perturbations off ( ρ, v i n i , ˜ v i , B i n i , ˜ B i , p, φ ) T respectively with ˜ v i , ˜ B i indicating the components of v i and B i orthogonal to the boundary.While the full characteristic decomposition is lengthy, the characteristic decomposition asso-ciated to B n and φ is particularly simple. In fact, the associated co-basis is given byΘ L ± = 12 B n (0 , , , , , , ± c − l ) , (16)i.e. involving essentially just B n and φ , hence the imposition of suitable boundary conditionswill be directly related to that of our discussion of the simplified model in section 2.1. Boundary conditions
Having the complete characteristic structure it is now rather straightforward to determinethe possible boundary conditions. A simple one is to set all incoming modes to zero but this,generically, will not be consistent with the constraint. To this end, as explained in section2.1, one must analyze the induced evolution for the constraint and obtain from it a recipe forwhat to provide to the incoming modes. In the general case discussed here, notice that Θ L ± isessentially the same that gives rise to equation (12) since is non-trivial only for the B n and λ components. As a result, one can follow the same strategy and formulate constraint preservingboundary conditions by enforcing˙ U = U L + (cid:18)
12 ( F φ − F n ) (cid:19) ≡ £ ( U ) , (17)where £ ( U ) denotes the (maximally dissipative) constraint preserving boundary conditiondefined by equation (17).
3. Numerical Tests
In this section we present the results of tests aimed to examine the solution’s behavior withthe different possible choices. To this end we constructed an implementation of the MHDequations (eqns. 13) on a 2-dimensional domain with Cartesian coordinates ( x, y ) ∈ [ − L, L ] × [ − L, L ]. To discretize the system we employ a set of techniques which guarantee the stability ofgeneric linear hyperbolic systems. These techniques are constructed to satisfy at the discretelevel, the conditions holding at the analytical one to prove well posedness of a problem through nergy estimates [24]. Thus we adopt discrete derivative operators satisfying summation byparts, a 3rd order Runge-Kutta time integration and add dissipation through a Kreiss-Oligertype operator consistent with summation by parts. Finally, in all the tests considered we adda homogeneous ‘atmosphere’ throughout the computational domain in the initial data (and a‘floor’ during the evolution) for the density and the pressure to prevent negative values fromarising in the simulation. We typically adjust the atmosphere (floor) values to be 10 − (10 − )times the maximum of the initial values of ρ and e respectively. Notice that this straightforwardimplementation is not expected to handle shocks or discontinuities, however, since our goal isto examine boundary conditions we restrict here to scenarios where these effects do not arisewithin the time of interest.For testing purposes we adopt a few different initial scenarios which are based on modifica-tions of two commonly employed initial data sets: Rotor test:
This is a variation of the MHD Rotor problem [28] which is smooth and allowsfor considering an additional flow field and initial constraint violation, ρ ( x, y ) =
10 if r ≤ r f ( r ) if r < r < r r ≤ r (18) v x ( x, y ) = ν x + − f ( r ) v o yr if r ≤ r − f ( r ) v o yr if r < r < r r ≤ r (19) v y ( x, y ) = ν y + f ( r ) v o xr if r ≤ r f ( r ) v o xr if r < r < r r ≤ r (20)and p = 1 (21) B x = B + B x (22) B y = 0 (23) φ = 0 (24)with r = 1 , r = 2, r = √ x + y , γ = 1 . v a constant and f ( r ) an interpolating function.For the general tests we adopt f ( r ) = ( r − r )( r − r ) . (25)However, for the convergence test we adopt a smooth interpolating polynomial and velocityprofiles as ( r ) = P ( r ) (26) v o = ˜ v o r / (1 + r ) (27)with P ( r ) a polynomial of order 5 such that P ( r ) = 1; P ( r ) = 0 and both first andsecond derivatives at r = { r , r } are zero. Thus, the corresponding initial data is at least C throughout the computational domain.Finally, the functions ( ν x , ν y ) and B x are included to consider further modifications for specifictests, in particular choosing :– B x = κe − r allows one to introduce data violating the constraints.– ( ν x , ν y ) = − ǫ (cid:16) − e − r (cid:17) ( x, y ) defines data inflowing to further test boundary issues. Blast test:
This is a variation of the MHD Blast problem [29].The initial data is given by, ρ = 1 . v x = 0 . v y = 0 . B x = 4 . B y = 0 . p = e − r / φ = 0 . U i = 0, for i = 1 .. + ˙ U = 0, withΞ + = P i U + j θ + j , for j such that λ j > + ˙ U = £ ( U ), as defined by eqn.(17) (denoted by CP).in all tests, and compare the solutions’ behavior when employing each option. We begin byconsidering the flux-conservative form of the equations (keeping α = 0) and then examineparticular cases with the Galilean invariant form ( α = 1). Testing the implementation
In the first test we confirm the overall convergent behavior of the numerical solution whenconsidering sufficiently smooth initial data. We evolve the C version of the Rotor initial data Time || ρ x2 - ρ x1 || || ρ x4 - ρ x2 || ||B x2 - B x1 || ||B x4 - B x2 || Fig. 1. Behavior of the L norm of the difference between solutions obtained with gradually improved resolutions. The differencebetween them converges to zero as expected. . (with no constraint violation or background fluid flow) for three different resolutions ∆ l =1 . / l ( l = 0 , ,
2) and check the pair-wise difference of the numerical solutions obtained de-creases as expected. Figure 1 illustrates the behavior of || F (∆ l ) − F (∆ l +1 ) || (with F either ρ or B x ). As is evident in the figure, as resolution is improved the differences decrease as ex-pected. For all remaining tests we present results for the finest grid employed with l = 2. Blast initial data
In this test, we adopt the blast initial data, evolve it for different choices of boundary condi-tions setting α = 0 , c l = 20 and s = 1 and examine the constraint’s behavior in each case. Fig-ure 2 shows the L norm of the constraint for the different boundary value conditions. Clearlythe numerical solution obtained with constraint preserving boundary conditions is superiorby about an order of magnitude in constraint violation than the no-incoming case and almostthree orders better than that obtained with the freezing boundary condition. An importantpoint to emphasize is that a closer inspection of the solutions obtained reveals that the maincontribution to the error originates at the boundaries in all cases (though with essentially thesame behavior as far as the error’s magnitude with respect to the boundary condition adopted).The norm displayed in Figure 2 is calculated over the whole computational domain ignoring the last two points at all boundaries to avoid placing excessive weight on the violation at theboundary. Nevertheless, as indicated, constraint preserving boundary conditions give rise to asolution whose constraint is violated the least. Rotor initial data. Effects of boundary conditions and divergence cleaning
We adopt this data to further examine the solution’s behavior and allowing for initial viola-tions of the constraint. To this end we adopt κ = 10 − . Time || D || NICPFR
Fig. 2. Behavior of the L norm of the constraint for the different boundary conditions. The violations in each case grow fromround-off values before settling to an approximately constant value. The values throughout the run for this norm improves asthe boundary condition adopted is refined as expected, in fact the solution obtained with the constraint preserving boundarycondition preserves the constraint by at least three-order of magnitude better than the one obtained with the freezing condition. . Figure 3 illustrates the solution’s behavior with the different boundary conditions and withand without the use of the divergence driver. The addition of the divergence driver allows fora dynamical reduction in the constraint violation until t ≃ .
6, at this time the propagatingmodes interact with the boundaries which become the main source of error. Here again onesees that the constraint preserving boundary condition gives rise to significantly smaller er-rors in the solution. Interestingly however, adopting the no-incoming modes together with theconstraint damping field provides a reasonably similar behavior. This is due to the fact thatthe no-incoming boundary condition allows the outgoing constraint violating mode to leavethe computational domain, while the incoming constraint violating mode, generated by notimposing the constraint preserving boundary condition, is damped to a significant degree bythe divergence cleaning technique.Another interesting behavior is revealed when varying c l . This affects the constraint violatingmodes’ propagation speeds which in turn has a strong influence in the solution’s constraintbehavior. In what follows we adopt the constraint preserving boundary condition. For thistest we use the Galilean invariant system, setting s = 0 (so as not to include damping) andchoose c l ranging from 10 to 80 ⋆ together with adjusting the time step accordingly in ordernot to violate the Courant condition. Figures 4 and 5 (with initial constraint violation) showthe effect of varying c l in two cases, the Rotor initial data and its modification to include abackground velocity given by ( ν x = ± . x, ν y = ± . y ). We concentrate on the latter casefor there the effect is more striking. Notice that, in figure 4, the initial plateau correspondsto round-off values since the constraint is initially satisfied to that level. Once the non-trivialpart of the solution reaches the boundary a significant violation of the constraint is generated. ⋆ For smaller values of c l ( c l <
10) the code became unstable due to instabilities generated at the corners of the computationaldomain. Time || D || CP φ FR φ INC φ CPFRNI
Fig. 3. Behavior of the L norm of the constraint for different options with some non-trivial initial violation of the constraints.As is evident in the figure, the divergence cleaning is able to damp the constraints by several orders of magnitude while theviolation is present in the bulk of the computational domain. After the solution reaches the boundary, the boundary valuesinduced there dominate the violation of the constraint and again the behavior is significantly improved by the no-incomingand constraint preserving conditions. Time || D || cl=20 cfl=0.025 wholecl=20 cfl=0.025 maskedcl=40, cfl/2 wholecl=40, cfl/2 maskedcl=80, cfl/4 wholecl=80, cfl/4 masked Fig. 4. Behavior of the L norm of the constraint for different values of the coupling constant c l with constraint preservingboundary condition. The norm is calculated over the whole computational domain (whole) or over an interior region of 3 / c l is increased a slight improvement is achieved. Depending on the boundary condition, and on the characteristic velocities of the system underconsideration that violation propagates to the inside or just leave the domain. In our case, theability of the constraint preserving boundary condition to allow constraint violating modesto leave the domain, results in smaller errors for faster propagating cases. This is a result ofviolations induced by the evolution leaved the computational domain more rapidly. Time || D || cl=20, cfl=0.025 wholecl=20, cfl=0.025 maskedcl=40, cfl/2 wholecl=40, cfl/2 maskedcl=80, cfl/4 wholecl=80, cfl/4 masked Fig. 5. Similar to figure 4 but with initial data containing constraint violations.
4. Conclusions
We have investigated the numerical solution to the (ideal) Newtonian magnetohydrodynam-ics equations and developed boundary conditions which at the continuum level are constraintpreserving and so would not introduce further violations in the computational domain. Atthe numerical level these conditions do introduce some truncation-error level violations whichconverge away with resolution. These boundary conditions developed are consistent with max-imally dissipative ones and so the discrete energy of the system remains bounded. We alsoexamined the constraint’s behavior when enlarging the system so as to couple an extra field.The addition of such field, together with a suitable modification of the equations induces anon-zero speed in the propagation of the constraint and allows for driving its value to zero.We have studied, in particular, two of the many equivalent systems. One which is fully conser-vative –which is only weakly hyperbolic without the addition of the extra field which rendersit symmetric hyperbolic– and a Galilean invariant one which is symmetric hyperbolic –evenwithout adding the extra field– but is not expressible in conservative form.In order to examine the solution’s behavior we implemented the equations in Cartesiancoordinates and with a spatial discretization that preserves the initial constraint error. Thisallows us to separate bulk from boundary effects and observe that the main violations do indeedoccur at boundaries. To examine the boundary condition effects on the evolution we adoptthree types of boundary conditions:– The most direct one, and crudest, sets to zero all right hand sides in a buffer boundary region.This might lead to inconsistencies for one might be prescribing conditions to outgoing modes.This condition is essentially equivalent the most commonly employed one of flux-copyingnear boundary points.– The second one is the “no-incoming” boundary condition where one projects out all incomingmodes in the evolution equations leaving the rest (tangential and outgoing modes) intact.The implementation of this condition is slightly delicate as modes can change directions overtime, and so may turn from being incoming to being outgoing and vice-versa. The third condition preserves the constraint in the sense that is deduced by setting to zerothe incoming mode related to the constraint propagation. This condition is a modificationto the previous one where now instead of projecting out all incoming modes, one incomingmode is fixed so that the incoming constraint mode is zero as a result.Not surprisingly the first option displays the worst behavior for the constraint as not onlynothing is done to minimize the constraint violations introduced through the boundaries butfurther inconsistencies in the solution are induced there. As a result, in the best case theconstraint mode bounces back from the boundary into the integration region. The secondoption performed substantially better displaying a gain of an order of magnitude in constraintpreservation. This is a result of the boundary condition allowing for one of the φ modes –carrying constraint violations– to leave the integration region. With the third alternative anadditional order of magnitude (at least) is gained with respect to the non-incoming condition.This condition not only allows for constraint violations to leave the computational domain butdoes not introduce significant violations at the boundary.Thus we see that appropriately handling the boundary conditions the problem of constraintbehavior can be significantly controlled. It should be mentioned that we still do not havea mathematically rigorous proof that the constraint preserving boundary conditions is wellposed. However, in simpler systems –like the toy model we discussed at the beginning of thepaper– a proof could be devised at least for the case where the eigenvalues do not change sign.Further work in this direction is needed to put the analysis and results obtained in strongergrounds.
5. Acknowledgments
We would like to thank M. Anderson, E. Hirschmann, D. Neilsen and C. Palenzuela, for manystimulating discussions during the course of this work. This work was supported by the NationalScience Foundation under grants PHY-0326311, PHY-0554793, 0653375 to Louisiana StateUniversity and PHY05-51164 to the Kavli Institute for Theoretical Physics. L.L. thanks theKavli Institute for Theoretical Physics for hospitality where parts of this work were completed.
Appendix A. Characteristic decomposition
In this appendix we revisit the characteristic analysis of the ideal MHD system coupled tothe divergence cleaning field (a related discussion in the absence of this field can be found in[30,31]). Using the decomposition on normal and tangential parts, and defining d := σ + v n weget: d ¯ ρ = − ρ ¯ v n ,d ¯ v n = 2 B n /ρ ¯ B n − /ρB j ¯ B j − /ρ ¯ p ,d ¯ B n = v n ¯ B n + ¯ φ ,d ¯ p = − c s ρ ¯ v n − ( γ − B j v j ¯ B n + (1 − γ ) B n ¯ φ ,dB ¯ v = B /ρ ¯ B n − B n /ρ ¯ p ,dB ¯ B = Bv ¯ B n − B ¯ v n + B n B ¯ v + B n ¯ φ ,d ˜¯ v i = ˜ B i /ρ ¯ B n + B n /ρ ˜¯ B i , ˜¯ B i = ˜ v i ¯ B n − ˜ B i ¯ v n + B n ˜¯ v i .σ ¯ φ = c l ¯ B n (A.1)For the sake of clarity we will present the solution to the eigenvalue/eigenvector problemalong different cases which are naturally divided by the behavior of certain modes. A.0.1.
BASIS • Normal modes . These correspond to ¯ B n = ¯ φ = 0, where the equations become: d ¯ ρ = − ρ ¯ v n ,d ¯ v n = − B ¯ B/ρ − ¯ p/ρ ,d ¯ p = − c s ρ ¯ v n ,dB ¯ v = − ( B n /ρ )¯ p ,dB ¯ B = − B ¯ v n + B n B ¯ v ,d ˜¯ v i = ( B n /ρ ) ˜¯ B i ,d ˜¯ B i = − ˜ B i ¯ v n + B n ˜¯ v i . (A.2)Setting d = 0 one can show that all but the first component must vanish while itself can haveany value. So a first eigenvector is given by U = (1 , , , , , , , (A.3)where the entries are: U = ( ¯ ρ, ¯ v n , ˜¯ v i , ¯ B n , ˜¯ B i , ¯ p, ¯ φ ).Using all the equations we get the eigenvalue condition: d − ( c s + B /ρ ) d + c s B n /ρ = 0 . (A.4)From which we get four solutions: d p ± = ± s(cid:18) ( c s + B /ρ ) + q ( c s + B /ρ ) − B n /ρ ) c s (cid:19) / d m ± = ± s(cid:18) ( c s + B /ρ ) − q ( c s + B /ρ ) − B n /ρ ) c s (cid:19) / B → d p solutions tend to ± c s , while the d m → k B n k / √ ρ . The d p solutions tend to the fluid solutions while the d m to pure magnetic ones, so they decouplein this limit. One thus must choose the eigenvectors’ normalization in a suitable way to reflectthis behavior. For the d p solutions we set ¯ v n = 1, and use the last two equations to computethe normal magnetic field and the normal velocity, obtaining: U P ± = ∓ ρd p , , − B n R p ρ ˜ B i , , ∓ d p R p ˜ B i , ∓ c s ρd p , ! , (A.7)where R p := d p − B n /ρ , has a finite limit as B → or the d m solutions we proceed the other way around, we fix B ¯ B = B i ˜¯ B i = B i M i , where M i = ˜ B i / | ˜ B | and so compute ¯ v n , and ¯ p from the second and third equations of system (A.2),obtaining: U M ± = B ¯ BR m , ∓ d m B ¯ BR m ρ , ± B n ρ d m M i , , M i , c s B ¯ BR m , ! , (A.8)where R m := d m − c s has a finite limit as B → B i = M i perpendicular to ˜ B i , but otherwise arbitrary,that is B ¯ B = ˜ B i ˜¯ B i = 0. In this case all other components vanish except the tangential velocitywhich can be expressed as ˜¯ v i = d a /B n A i with ~A = ( − ˜ B , ˜ B ) / | ˜ B | and d a ± = ± B n / √ ρ (sothat the last two equations can have a nontrivial solution). The corresponding eigenvector is: U A ± = , , ± d a B n A i , , A i , , ! . (A.9) • φ modes We now look at the modes which are induced by the introduction of the field φ . There, fromthe third and last equation of (A.1) we get a 2x2 system which implies: σ = ± c l and so, d := d l = ± c l − v n . If we fix ¯ B n = b we can obtain the scalar components by solving the followingsubsystem of (A.1): d l ¯ v n + B ¯ B/ρ + ¯ p/ρ = 2( B n /ρ ) b ,c s ρ ¯ v n + ¯ pd l = ( − ( γ − Bv + (1 − γ ) σB n ) b , − B ¯ v n − d l B ¯ B + B n B ¯ v = − ( B n σ + Bv ) b , ( B n /ρ )¯ p + d l B ¯ v = ( B /ρ ) b . (A.10)From this subsystem we get:¯ v n = [ − ( d l − B n /ρ ) F + d l F + d l ρ ( d l F − ( B n /ρ ) F )] /δ l , (A.11)¯ p = ( − c s ρ ¯ v n + F ) /d l , (A.12)where: δ l = ρ ( d l − d l ( c s + B /ρ ) + c s B n /ρ ) ,F = − B n bρ ,F = (1 − γ )( Bv + B n σ ) b ,F = − F / (1 − γ ) ,F = B bρ . Once we have ¯ v n we can solve for the vectorial components of the system: l ˜¯ v i − ( B n /ρ ) ˜¯ B i = ˜ B i b/ρ , (A.13) d l ˜¯ B i − B n ˜¯ v i = ˜ v i b − ˜ B i ¯ v n , (A.14)obtaining:˜¯ v i = [˜ v i bB n /ρ + ˜ B i ( d l b − ¯ v n B n ) /ρ ] /δ s , (A.15)˜¯ B i = [˜ v i d l b + ˜ B i ( − d l ¯ v n , + bB n /ρ )] /δ s , (A.16)where δ s := d l − B n /ρ . Combining the above intermediate steps we finally obtain, U L ± = − ( ρ/d l )¯ v n , [ − ( d l − B n /ρ ) F + d l F + d l ρ ( d l F − B n F ρ )] /δ l , [˜ v i bB n /ρ + ˜ B i ( d l b − ¯ v n B n ) /ρ ] /δ s , b, [˜ v i d l b + ˜ B i ( − d l ¯ v n + bB n ρ )] /δ s , ( − c s ρ ¯ v n + F d l , bc l /σ ! (A.17) A.0.2.
CO-BASIS
Net we compute the co-basis which we will use to construct the suitable projector for en-forcing different boundary conditions. As in the previous analysis, we divide our task alongdifferent eigenvalues for simplicity. • Θ Since U has only one non-vanishing component (the first one) all elements of the co-base,except the corresponding one have a vanishing first element.The first element, is simply:Θ = (1 , , , A, , − c − s , B ) . (A.18)To find the remaining elements we notice that if we define: V P := (
U L + + U L − ) and V M :=(
U L + − U L − ), the first has a zero in the last component (corresponding to ¯ φ ) while the secondhas a zero in the fourth component (corresponding to ¯ B n ). Thus contraction of Θ with V P and
V M will leave either A or B only as unknowns. Thus we might proceed to setting A and B temporarily to zero, perform the contraction and extract what these must be obtaining, A = − Θ ( V P )2 b , (A.19) B = − Θ ( V M )2 bc l , (A.20)where we first temporarily set to zero A , and B in Θ . • Θ P ± From P ± ( U M + − U M − ) = 0 andΘ P ± ( U M + + U M − ) = 0we get that the following structure for it:Θ P ± = , C, C d m ˜ B i B n R m , D ± , − E ± c s ˜ B i R m , E ± , F ± ! . (A.21)Now, from 1 = Θ P ± ( U P + + U P − ) we obtain: C = R p R m R p R m − d m ˜ B /ρ ) (The same for both). (A.22)While from ± P ± ( U P + − U P − ) we deduce E ± = ± R p R m d p c s ( d p ˜ B − ρ R p R m ) . (A.23)The remaining two components can be computed using V P and
V M as above obtaining, D ± = − Θ P ± ( V P )2 b , (A.24) F ± = − Θ P ± ( V M )2 bc l , (A.25)where we set to zero D , and F in Θ P ± , as in the previous case. • Θ M ± From Θ M ± ( U P + − U P − ) = 0 and Θ M ± ( U P + + U P − ) = 0 we get that the following structurefor it:Θ M ± = , B n ˜ B i N i ± ρR p , N i ± , S ± , L i , − d p ˜ B i L i c s R p ρ , T ± . (A.26)From 1 = Θ M ± ( U M + + U M − ) we obtain: L i = − R m R p M i
2( ˜ B d p /ρ − R m R p ) . (A.27)From ± M ± ( U M + − U M − ) we obtain: N i ± = ∓ d m R m R p M i B n /ρ ( ˜ B d m /ρ − R m R p ) . (A.28)The remaining two components, S and T are computed as in the previous case, using: V P and
V M removing the ¯ B n and the last component ¯ φ in Θ M ± . We get, S ± = − Θ M ± ( V P )2 b , (A.29) T ± = − Θ M ± ( V M )2 bc l . (A.30) Θ A ± We get the co-vector of
U A by direct computation:Θ A ± = , , ± B n A i d a , G ± , A i , , H ± ! , (A.31)where the components, G and H are given by: G ± = − Θ A ± ( V P )2 b , (A.32) H ± = − Θ A ± ( V M )2 bc l . (A.33)removing the components ¯ B n and ¯ φ in Θ A ± . • Θ L ± It remains now to determine the last two elements. Since the first 7 eigenvectors span completelythe seven dimensional space given by ( ¯ B n = ¯ φ = 0) the co-vectors have only components inthat subspace and are given by:Θ L ± = (cid:18) , , , b , , , ± bc l (cid:19) . (A.34) Caveat: Singular points
The characteristic structure outlined above may change as some eigenvalues can change mul-tiplicity for particular values of the fields. These special cases require further analysis. Recallthe eigenvalues are given by: σ = − v n ,σ + P = d p − v n ,σ − P = − d p − v n ,σ + M = d m − v n ,σ − M = − d m − v n ,σ + A = d a − v n ,σ − A = − d a − v n ,σ + L = c l ,σ − L = − c l , they can cross when d p = 0, v n = ± c l , d m = 0 and d p = d m . The first two cases can be dealtwith by appropriately chosen parameters so that they would not occur in physically relevantscenarios. For instance, the first case would imply a vanishing sound speed and will not beconsidered since we will not deal with a fluid describing dust. The second case could be avoidedby simply taking c l large enough. Therefore, we concentrate on the other two: • d m = 0 or this we need that c s + B ρ ! = c s + B ρ ! − c s B n ρ (A.35)which only happens when B n = 0. Near that point we have, d m = c s B n ρc s + B . (A.36)At this point we compute explicitly the eigenspace for the coinciding eigenvalues and chooseeigenvectors with good limit. We use them in a sufficiently small neighborhood of this pointusing conditional statements in the code at every point in the boundary.In three dimensions, in this case, also d a = 0, therefore the eigenvectors U M ± and U A ± aredegenerates. • d m = d p For this to happen we need0 = ( d p − d m ) = c s + B ρ ! − c s B n = c s − B n ρ ! + 2 c s ˜ B i ρ + ˜ B i ρ . (A.37)Thus. B n /ρ = c s , and ˜ B i = 0. Notice that in that case we have R p = R m = 0. We proceed in asimilar fashion as above. Notice that when eigenvalues coincide all what enters in the boundarycondition is the projector on that subspace, so one just has to choose the eigenvectors so thatnumerically that projector is robust in a small neighborhood.In the case of three dimensions, d a = B n /ρ . Hence, d a = d m = d p and the eigenvectors U P ± , U M ± and U A ± become degenerates.References [1] H. Goedbloed and S. Poedts, “Principles of Magnetohydrodynamics”. Cambridge Univ. Press, Cambridge, 2004.[2] P.A. Davidson, “An introduction to Magnetohydrodynamics”. Cambridge Univ. Press, Cambridge, 2001.[3] J. F. Hawley and C. Evans, “Simulation of mhd flows: A constrained transport method,” Astrophys. , vol. 332, p. 659, 1988.[4] G. Toth, “The ∇ B constraint in Schock-Capturing Magnetohydrodynamics”, J. Comp. Phys. , 605-652 (2004).[5] C. F. Gammie, J. C. McKinney and G. Toth, Astrophys. J. , 444 (2003)[6] T. A. Gardiner and J. M. Stone, “An unsplit godunov method for ideal mhd via constrained transport,” J. Comput. Phys. ,vol. 205, pp. 509–539, 2005.[7] S. A. Balbus and J. F. Hawley, “Numerical simulations of mhd turbulence in accretion disks,”
Lect. Notes Phys. , vol. 614,pp. 329–350, 2003.[8] M. Shibata and Y.-i. Sekiguchi, “Magnetohydrodynamics in full general relativity: Formulation and tests,”
Phys. Rev. ,vol. D72, p. 044014, 2005.[9] M. D. Duez, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, “Relativistic magnetohydrodynamics in dynamical spacetimes:Numerical methods and tests,”
Phys. Rev. , vol. D72, p. 024028, 2005.[10] B. Giacomazzo and L. Rezzolla, “Whiskymhd: a new numerical code for general relativistic magnetohydrodynamics,”
Class. Quant. Grav. , vol. 24, pp. S235–S258, 2007.[11] S. Li, H. Li and R. Cen, arXiv:astro-ph/0611863.[12] M. Obergaulinger, M. A. Aloy and E. Muller, Astron. Astrophys. , 1107 (2006)[13] J.S. Kim, D. Ryu, T.W. Jones and S.S. Hong, Astrophysical J. , 506 (1999).
14] A Dedner, F Kemm, D Kr¨oner, C-D Munz, T Schnitzer, and M Wesenberg. Hyperbolic divergence cleaning for the MHDequations. J. Comput. Phys., 175:645, 2002.[15] D.S. Balsara and J.S. Kim, Astrophys. J. , 1079 (2004).[16] M. Anderson, E. Hirschmann, S. L. Liebling and D. Neilsen, Class. Quant. Grav. , 6503 (2006).[17] D. Neilsen, E. W. Hirschmann and R. S. Millward, Class. Quant. Grav. , S505 (2006)[18] A. Dedner, D. Kroner, I. Sofronov and W. Wesenberg, J. Comput. Phys. , 448 (2001).[19] R. Leveque, “Finite Volume Methods for Hyperbolic Problems”, Cambridge University Press, 2002.[20] A. Zachary, A. Malagoli and P. Colella, SCIAM J. Sci. Comput, , 263 (1994).[21] D. Balsara, Astrop. J. Supp. , 119 (1998).[22] D. De Zeeuw, T. Gombosi, C. Groth, K. Powerl and Q. Stout, IEEE Trans. Plasma Science, , 1956 (2000).[23] H. Yee and B. Sjogreen, in Proceedins of the International Conference on High Performance Scientific Computing, Hanoi,Vietnam. (2003). Cambridge University Press, 2002.[24] B. Gustaffson, H. Kreiss, and J. Oliger, Time Dependent Problems and Difference Methods (Wiley, New York, 1995).[25] Oscar A Reula, Journal of Hyperbolic Differential Equations 1:22, 251-269, 6/2004[26] S. K. Godunov, The symmetric form of magnetohydrodynamics equation, Numer. Methods Mech. Contin. Media 1, 26(1972).[27] T. J. Barth, Numerical methods for gasdynamic systems on unstructured meshes, in An Introduction toRecent Developments in Theory and Numerics for Conservation Laws: Proceedings of the International School,Freiburg/Littenweiler, Germany, October 20?24, 1997, edited by D. Kr¨oner, M. Ohlberger, and C. Rohde, Lecture Notesin Computational Science and Engineering (Springer-Verlag, Berlin, 1999), Vol. 5, p. 195.[28] D. Balsara and D. Spicer, Journal of Comp. Phys., , 397, (2003).[30] S. Pennisi, “A covariant and extended model for relativistic magnetofluiddynamics”. Ann. Inst. Henri Poincare , 343-361(1993).[31] D. Balsara, “Total variation diminishing for relativistic magnetohydrodynamics”. The Astrophysical Journal Supp. Series, , 83-101 (2001)., 83-101 (2001).