ture MHD shocks and contact discontinuities, and to model accurately Lorentz force ... 1National Center for Supercomputing Applications, University of Illinois at ...
A Test Suite for Magnetohydrodynamical Simulations James M. Stone1;2, John F. Hawley3, Charles R. Evans4;5 and Michael L. Norman1;2 ;
Received:
Accepted:
Abstract We have drawn together in this paper a collection of magnetohydrodynamic (MHD) problems that will provide researchers who have interests in modeling MHD phenomena with a battery of tests (a \test suite") for calibrating their numerical algorithms. Use of these tests will provide a common reference for comparison of di erent numerical MHD algorithms. The test suite includes both one and two dimensional problems. Taken together, these problems test the abilities of a numerical method to propagate accurately all of the MHD wave families in moving and stationary media, to capture MHD shocks and contact discontinuities, and to model accurately Lorentz force terms in multiple dimensions. An example solution of each test problem is presented. We demonstrate diagnostic convergence-testing procedures that provide a quantitative evaluation of MHD algorithms.
Subject headings:
1
hydrodynamics|MHD
Introduction
Although it has long been recognized that magnetic elds play important roles in at least a few astrophysical environments, a growing body of observational evidence suggests 1 National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign 2 Also Department of Astronomy, University of Illinois at Urbana-Champaign 3 Department of Astronomy, University of Virginia 4 Department of Physics and Astronomy, University of North Carolina 5 Presidential Young Investigator; Alfred P. Sloan Foundation Fellow AAS WGAS macros v2.0
1
that magnetic elds are both ubiquitous ( Beck 1990) and may in many instances dominate the dynamics of astrophysical systems. For example, magnetic elds probably control the evolution of the interstellar medium, particularly in dense molecular clouds and star-forming regions ( Mouschovias 1987; Shu, Adams, & Lizano 1987). Magnetic elds are thought to be important in determining the dynamics of gas ows near the Galactic Center (Wardle & Konigl 1990), in accretion ows and disks, and in extragalactic radio jets ( the review by Blandford 1989 and the proceedings edited by Belvedere 1989). Detailed magnetohydrodynamic (MHD) models are needed to investigate these many scenarios. However, the complexity of the MHD equations guarantees that most detailed modeling will rely on numerical simulations. e.g.
et al
e.g.
e.g.
Codes designed for large-scale MHD simulations must be thoroughly tested to eliminate errors in the implementation of a given algorithm and to check the accuracy of the algorithm. Usually, a code test consists of comparing a numerical solution to either an analytic solution or to a numerical solution computed with another code of established performance. Typically, analytic solutions are so simpli ed that rarely does one problem test all of the terms in the full MHD equations. A set of problems is required to test a code thoroughly and we refer to such a collection as a \test suite." Recently, during the course of investigating numerical MHD algorithms as part of development e orts on the general purpose radiation MHD code ZEUS-2D (described in Stone & Norman 1991a,b; and Stone, Mihalas, & Norman 1991) and a separate MHD and gravity code (Evans & Hawley 1988), we identi ed and subsequently used a comprehensive set of MHD test problems. Because of their potential importance to other researchers, we have assembled this set of tests in this paper. Each problem in the test suite focuses on a di erent physical e ect, including uid and magnetic eld advection, propagation of linear and nonlinear MHD waves, and force balance in multiple dimensions when the Lorentz force, J 2 B , is dynamically signi cant. 2
Individual tests and sets of tests for hydrodynamic codes have been presented before. For example, Sod (1978) introduced a Riemann shock tube problem that has become a standard test for hydrodynamic codes. Unlike Sod (1978), this paper is not primarily intended to report the results of the calibration of the algorithm used in any particular code (although we present a representative numerical solution for each test problem computed with our method). Nor do we compare a wide variety of MHD algorithms; numerical MHD is insuciently developed at present to make this possible. Instead, this paper de nes and describes the test problems, and provides a reference numerical solution for each, with the anticipation that they will nd wide use as additional MHD algorithms are developed. A second major emphasis of this paper is the use of more sensitive test diagnostics. As the use of numerical techniques becomes more widespread in astrophysics, algorithm tests must become quantitative; it is no longer adequate merely to plot the numerical result over the analytic solution and quote \representative" percentage errors ( Hawley, Smarr & Wilson 1984). It is far more useful and important to plot errors directly and compute convergence rates as a function of mesh resolution. The simple prescriptions given here provide a detailed quantitative assessment of code and algorithm performance, thus signi cantly increasing the ecacy of the testing process. Hence besides providing a standard set of test problems, we anticipate that this test suite will provide a true quantitative baseline for comparing other numerical algorithms. e.g.
The organization of this paper is as follows. We begin in x2 by discussing the equations of ideal MHD, the important aspects that any numerical algorithm should incorporate in solving these equations, and the most e ective methods by which test problems can be used to calibrate these algorithms. In x3 we de ne each of the ve test problems in our test suite, and give representative numerical solutions. Finally, in x4, we present our conclusions.
3
2
Numerical Magnetohydrodynamics
The equations of ideal MHD (in cgs units) are D + r1 v = 0; Dt
Dv Dt
= 0rp 0 r8 + 41 (r2 B) 2 B; !
D e = 0pr1 v; Dt @B = r2 (v 2 B ); @t r28 = 4G:
(1) (2) (3) (4) (5)
where the dependent variables are the mass density , the velocity v, the internal energy density e, the magnetic eld B , and the gravitational potential 8. The D=Dt denotes the Lagrangian derivative, and the system of equations is closed by specifying an equation of state p = p(; e). The test problems described below all assume ideal MHD so that the magnetic eld is assumed to be \frozen-in" to the uid. Of course, in some applications ideal MHD will not apply. A good example is dense molecular clouds where the degree of ionization of the plasma is low enough that neutral gas drifts relative to charged particles (ambipolar di usion). To model systems exhibiting ambipolar di usion, or any other transport e ect not accounted for in equations [1]-[5] ( resistivity), additional terms must be added, and additional problems speci cally designed to test these new e ects would be required. e.g.
The numerical solutions of equations [1]-[5] that are presented in this paper use an algorithm described in Stone and Norman (1991b), to which the interested reader is referred for details. Here we merely point out the guiding principles behind this scheme. The algorithm is based upon the constrained transport (CT) formalism (Evans & Hawley 1988) for directly evolving the magnetic eld while maintaining the divergence-free constraint. The CT formalism is the nite-di erence implementation of Stoke's theorem; it relies upon the computation and consistent use of electromotive forces (EMFs). The CT scheme places no 4
particular requirements on how the EMFs are to be computed. In fact, therein lies the second cornerstone of our MHD method: we use the method of characteristics (MOC) to compute both the EMFs used in the induction equation [4] and the transverse component of the J 2 B force in equation [2]. The MOC-CT procedure propagates Alfven waves with less dispersion, allowing stable treatment of discontinuous waves, abrupt changes in propagation velocity, and outgoing and re ecting boundary conditions. More obvious di erencing techniques were shown to yield unacceptable results. Useful alternatives to MOC-CT may exist and the test problems described in this paper provide a means to experiment with such algorithms. An important facet in testing any algorithm is the comparison between the numerical and reference solutions. Typically numerical results are plotted on top, or along side, of the reference solution. At best, however, this technique only allows a qualitative check on the numerical results; there are several distinct disadvantages to this practice. For one, the presence of the \analytic" curve in the plot can lead to the perception that the numerical solution contains more structure than is really there. A good example is the blast wave solution in spherical geometry, which features a sharp cusp-like maximum in density and pressure immediately behind the shock, with density and pressure rapidly decreasing in the post-shock region. Typical numerical solutions smear the shock out over several zones and cannot resolve or capture the cusp. The maximum numerical value is typically substantially less than the analytic, yet the overlaid analytic plot can incorrectly lend the impression that the ne structure has been adequately treated. A second drawback involves the typical small size of errors and the attendant diculty in judging di erences even qualitatively. The solution to both issues is to plot the errors themselves. With attention directed to the errors, the next step becomes apparent: evaluating how the errors change with numerical resolution. In particular, the rate at which the error decreases with resolution (the ) can be measured. The issue of convergence is a central consideration in formal numerical analysis but it can also be applied on an experimental basis in re ning numerical algorithms (see for example Finn & Hawley 1989 convergence rate
5
for application of convergence testing to numerical advection). In addition, the convergence properties of numerical solutions relative to a reference is a powerful tool for detecting subtle coding errors. In developing the MHD test suite we have made use of three techniques for measuring the convergence rate. These are (1) convergence to an analytic solution measured at a xed time, (2) self-convergence measured at a xed time, and (3) single timestep convergence testing. We describe each of these in turn. We de ne the relative error at some zone labeled by i as i = jqi 0 q~(xi )j=q~(xi);
(6)
where qi is the numerical solution and q~(xi) the reference solution. The absolute error is similarly de ned except that the di erence is not normalized by q~(xi). We obtain a global measure of the error by taking a norm. For example, the L1 error norm is N X 1 i ; (7) L1 = N
i=1
where N is the number of gridpoints. This error should be asymptotically proportional to some power of the discretization interval, L1 4xr , where r is the convergence rate, as 1x ! 0 for a xed ratio 1t=1x. This corresponds to the formal order of the scheme. At any nite resolution, however, the error may converge faster or slower than the asymptotic value. If an analytic solution is available, the rst of the above mentioned convergence tests can be exploited, and the convergence of an algorithm to the exact solution measured. Analytic initial conditions are placed on the grid and evolved to a xed time. The global error is computed from the analytic solution to obtain an error norm. This procedure is repeated for various resolutions and the change in the value of the error norm as a function of resolution gives the convergence rate (slope of the log 0 log graph). Note that the rst thing to check, however, is whether or not the numerical solutions converge at all. There are many possible reasons for a failure to converge to an analytic solution, the most obvious of 6
which is the presence of a code \bug." In this case convergence testing is quite sensitive for it uncovers bugs in \unimportant" terms, i.e., those that have no apparent in uence on the solution values (but not the solution ). Convergence will not occur unless all of the nite-di erence equations are consistent with the analytic solution. A primary example is if the numerical implementation of boundary conditions is not consistent with the analytic solution. errors
Assuming that the numerical solution can be made to converge to the analytic result, the convergence rate can be measured. It should be noted that without explicit measurement of convergence rate, claims of \nth-order accuracy" for a given code or algorithm are of little value. In the absence of an analytic solution, self-convergence testing can be used to analyze a code. However, caution is in order, since a code will generally converge to some solution even when there are errors in the implementation, although the solution is not the one intended. In any case one can still test convergence by computing a problem out to some xed time using several di erent resolutions and measuring the relative changes in all variables. A convergence rate can be computed by using at least three di erent resolutions and determining the change in the errors. It should be noted that long evolutionary times are not required for complete code testing. For example, one can carry out single timestep convergence testing with an analytic solution. The known solution is initialized on the grid and evolved for one timestep. An error norm for the evolved solution is computed and used to obtain an error per unit time. This is a particularly powerful technique since a large range of resolutions can be tested with a minimal investment of computer time. This can be especially useful in testing multidimensional codes. In addition, coding errors are more easily spotted since in one timestep a single implementation error will generally only appear as nonconvergence of one variable. By plotting errors and convergence rates, a numerical algorithm can be debugged much 7
more eciently. For this reason, we present both a representative numerical solution for each test problem, and whenever possible the errors and convergence rate of our numerical algorithm. The power of these testing techniques can in principle be even more fully exploited and the tests presented here may suggest the development of additional diagnostic procedures. 3
The Test Problems
In this section we describe the ve problems comprising the MHD test suite. They are (1) Eulerian advection, (2) an MHD Riemann problem, (3) the propagation of noncompressive transverse MHD waves (Alfven waves), (4) a 1-1/2 dimensional evolution of a stationary ow, and (5) a fully 2-1/2 dimensional dynamical problem with shocks. These problems represent a subset of the tests that we used in developing our methods and exhibit what we consider to be the most important features. Each problem tests a di erent MHD phenomenon. The rst tests the advection of a magnetic uid, the second tests the ability of the scheme to propagate nonlinear compressive waves (shocks) and contact discontinuities, the third tests the propagation of non-compressive transverse waves, and the last two are dynamical problems that test the overall implementation and performance of the scheme. While no number of code tests can guarantee code performance in all possible applications, this suite represents the minimal set of essential tests. We have found each test to be equally important; a scheme may yield an acceptable result for one test yet fail completely on another. 3.1
An Advection Test
An important characteristic of any Eulerian uid dynamics algorithm is its ability to advect distributions of the dependent variables with as little numerical dispersion and di usion as possible. This rst problem is designed to verify the ability of the MHD algorithm to advect distributions of magnetic eld. 8
Typically, advection tests are performed in one dimension with the advection velocity assumed constant and the dependent variable distributed in either a discontinuous or continuous fashion. We begin with the one-dimensional advection of a square pulse of transverse magnetic eld in Cartesian geometry, identical to the test found in Evans and Hawley (1988). However, here we use the tools of convergence testing to explore more thoroughly the rami cations of this simple problem. This test measures the ability of an algorithm to advect contact discontinuities and hence measures the numerical dispersion and di usion of the scheme. The current density, 4J = r2 B, is examined to determine whether the scheme introduces anomalous extrema or sign reversals in the current. This is an important consideration for MHD calculations since the current density couples directly to the dynamics through the Lorentz force. We consider a pulse that is fty zones wide, take the velocity to be constant and advect the pulse for a distance ve times its width. A Courant number of 0.5 is used for this test. For illustration, we compare the advection of magnetic eld using three di erent monotonic interpolation methods: piecewise-constant \donor cell," piecewise-linear with van Leer's monotonicity constraint (van Leer 1977), and the piecewise-parabolic advection (PPA) scheme, which is a component of the Piecewise-Parabolic Method for hydrodynamics (Colella & Woodward 1984). These three advection schemes have previously been tested in the ZEUS-2D code (Stone & Norman 1991a). The results of this advection test are shown in Figure 1. The donor cell scheme is the most di usive of the three schemes with discontinuities being spread over about fty zones. The van Leer scheme is much less di usive and discontinuities are held to about fourteen zones. The PPA scheme is the least di usive, producing discontinuities that are seven zones wide. We have also considered PPA with a contact-discontinuity steepener (Colella & Woodward 1984), which produces discontinuities only two zones wide. In Figure 2 we plot the associated current density for each scheme. Analytically, the cur9
rent density is a delta function located at the pulse edges, but the di usivity of each scheme spreads the current over a nite width with a nite amplitude. Because the current density, via the Lorentz force, has an important e ect on the dynamics, the advection algorithm must be able to maintain its pro le accurately. In our tests we nd that the donor cell scheme spreads the current density over the whole pulse width with low amplitude. In contrast, the PPA scheme produces a current density spike only seven zones wide with a maximum amplitude that is seven times larger than the donor cell result. The van Leer scheme produces an intermediate result. We also note that for all three schemes, the current density is monotonic with no anomalous sign reversals. As discussed in Evans and Hawley (1988), this is because the CT algorithm evolves the magnetic eld components directly. If instead a scheme is employed where the vector potential is evolved (to maintain the divergence-free constraint), then a second derivative must be evaluated to compute the current. As Evans and Hawley (1988) point out the monotonicity constraints, used in the van Leer and PPA schemes, can produce sign reversals in numerically computed second derivatives. Figures 1 and 2 provide only a qualitative measure of the accuracy of the three advection schemes. Quantitative comparison is made in Figure 3, where the relative error of each numerical solution is plotted. For all three schemes, the largest error occurs at the pulse edges, as expected. However, the PPA scheme clearly gives the most accurate result, both in terms of the amplitude and spatial extent of the errors. Computing the convergence rate of the L1 error norm de ned by equation [6] is only partially useful for problems involving discontinuities. One can establish whether a scheme is converging at all, but otherwise all schemes will converge at a rst order rate on problems of this type. This occurs because the number of zones used by each scheme to represent a discontinuity is independent of resolution, thus the error in the transition zones remains constant. The spatial thickness of the discontinuity decreases as 4x so the L1 error norm simply converges at rst order for all schemes as N increases. However, the overall amplitude of the error norm remains a signi cant point of comparison between di erent schemes. 10
We can reproduce the formal convergence rate of each advection scheme by considering problems which do not contain discontinuities. Examples of such tests are given by Stone and Norman (1991a), and Finn and Hawley (1989) who compute error norms for a variety of advection schemes using several test problems. To summarize these results, donor cell, van Leer, and PPA converge at rst, second and third order rates when tested on problems involving constant-velocity advection of smooth, continuous data in one dimension and with constant zone size 4x. In general one will not obtain these formal convergence rates in dynamic multi-dimensional simulations. 3.2
An MHD Riemann Problem
The next problem in the test suite is a one-dimensional, fully dynamical problem involving the propagation of nonlinear compressive waves. A frequently used test of this type in hydrodynamics is the shock tube of Sod (1978), against which many numerical algorithms have been compared. A good candidate for the MHD analogue of the Sod shock tube is the magnetic Riemann problem introduced by Brio and Wu (1988; hereafter BW) to demonstrate the non-convexity of the equations of MHD. For our purposes it is sucient to note that the non-convexity of the MHD equations implies that di erent modes of the same MHD wave family can propagate with identical velocities. Thus, in an MHD ow it is possible to nd compound waves that consist of a shock wave attached to a rarefaction fan of the same family ( a slow shock attached to a slow rarefaction fan, ). This situation does not occur for waves in hydrodynamics. The MHD Riemann problem therefore tests the ability of the numerical scheme to represent shocks, rarefaction fans, and compound waves in MHD
ows. e.g.
etc.
Like the Sod shock tube, the MHD Riemann problem is initialized at t = 0 with two uniform states separated by a discontinuity. BW made speci c choices for the left and right states. For a = 2 gas, the pressure, density, and x and y components of the magnetic eld in the left state are chosen to be pl = 1:0, l = 1:0, (B x)l = 0:75 and (B y )l = 1:0 11
respectively, while in the right state they are chosen to be pr = 0:1, r = 0:125, (B x)r = 0:75, (B y )r = 01:0. Except for B x, which is always constant, these initial values are summarized in Table 1 in the rst and last columns. The problem is initialized on a domain of x 2 (0; 100) using 800 zones, with the discontinuity initially located at the zone interface at x = 50. The problem is run to time t = 10. Note that the choices for the dynamical variables are identical to those in the Sod shock tube problem, except for the ideal gas constant and the presence of a magnetic eld with a kink at the discontinuity. The angle that the magnetic eld lines make with the shock normal, in this case given by tan = 4=3, is an important parameter. The solution to the MHD Riemann problem will be entirely di erent for initial elds that are parallel to the shock front, perpendicular to it, or at arbitrary angles between these extremes. Unlike the Sod shock tube, the MHD Riemann problem possesses no known analytic solution. BW present numerical solutions to this problem using four methods: a LaxFriedrichs scheme, a Lax-Wendro scheme, a ux-corrected transport scheme, and a second order upwind Roe-type scheme designed for this problem. BW describe the implementation of the last method in detail; we consider it to give the best available reference solution for this problem (as provided in gure 2 of BW). Figure 4 shows the MHD shock tube computed with the MOC-CT scheme; pictured are the density, pressure, x0 and y0components of velocity, and the y0component of the magnetic eld. The solution consists of, from left to right, a fast rarefaction fan (FR), a slow compound wave (SC), a contact discontinuity with values on the left (CDl) and right (CDr ), a slow shock (SS), and another fast rarefaction fan (FR). Graphically, our results compare very favorably to the results given by BW, and also by a more complex Godunov scheme (Zachary & Colella 1989). The quantitative values obtained in the run are listed in Table 1 together with all such values given by BW for comparison. One of the goals of this paper is to get beyond the qualitative comparison a orded by plots such as gure 4. Although no analytic solution exists for the MHD Riemann problem, one can still measure the errors and convergence rate of the numerical solutions by using 12
self-convergence testing. To determine explicitly the order of the truncation error for the scheme, we employ a uniformly spaced nite di erence grid with resolution 1x = h. Denote the result of the numerical evolution of the initial data from time t = 0 to time t0 on the nite di erence grid by h (x; t0). Assume that the nite di erence scheme converges to the exact solution 0(x; t0) as h ! 0. De ne (h; x; t0) h(x; t0) 0 0(x; t0) as the truncation error (more precisely the discretization error). For suciently small h, let the truncation error be proportional to hk . Thus h (x; t0) 0 h (x; t0) ' ( k 0 1)(h; x; t0):
(8)
We can now determine the order of the scheme by examining h at three or more di erent values of h. To do this for the BW shock tube, we evolve a set of runs to a xed time using in each case uniform time steps that are proportional to the spatial zoning (the spatial resolutions are powers of two between 64 to 1024 zones). We calculate an L1 error norm de ned as the di erence between the solution at resolutions h and 2h: " R 2h h # j R0 j dx : (9) L 1
dx
The results of this test are given in Table 2 and plotted in Figure 5. As can be seen from the table, the self-convergence rate of the L1 error norm for density is 1.0. This holds for the other variables as well. Since this problem contains discontinuities, this convergence rate is anticipated for the reasons described in x3.1. In Figure 6, we plot the relative errors in the density as computed at resolutions of 200 and 800 grid zones. The relative error is computed with equation [6] and by using a 1600 grid point computation as the reference solution. (Note that in the absence of an analytic solution, convergence must be measured against a given code's own high resolution simulation.) In the MHD Riemann problem the errors are largest at the shocks, the contact discontinuity, and the rarefaction fans. A comparison of plots such as Figure 6 with similar 13
plots obtained with other MHD algorithms would provide quantitative information about the relative performance of the algorithms. 3.3
Magnetic Braking of an Aligned Rotator
MHD ows also support non-compressive (transverse) Alfven waves. As these are unique to MHD, it is essential to have problems that verify the ability of the numerical scheme to propagate Alfven waves. The analytic solution for the magnetic braking of an aligned rotator given by Mouschovias and Paleologou (1980, hereafter MP) is well suited to this purpose. Physically, the problem considered by MP is the time evolution of a rotating cylinder of dense gas embedded in a homogeneous ambient medium and threaded by a uniform axial magnetic eld. The shear at the disk surface generates Alfven waves which propagate along the magnetic eld lines into the ambient medium (accelerating it) and into the disk (decelerating it). MP give an analytic description of the time evolution of the rotational velocity and toroidal component of the magnetic eld in both the disk and ambient medium. Although the original application of the solution given in MP involved torsional waves in cylindrical geometry, the propagation of shear Alfven waves in Cartesian geometry is identical and the same analytic solution applies. By Galilean invariance, the analytic solution is also applicable in a moving background. To initialize the test, we follow MP and choose the unit of length to be the half thickness of the disk (Z ) and the unit of time to be the Alfven crossing time in the external medium Z1x=2=B , where x is the density in the external medium, and B is the axial magnetic eld strength. Convenient choices are then Z = x = B = 1. A single dimensionless parameter, = d =x , where d is the density in the disk, determines the magnetic braking time of the disk and characterizes the solutions. For = 1, the impedances of the disk and external medium are matched and the solution consists of two pulses in the toroidal components of velocity and magnetic eld which propagate out of the disk and to in nity in the external 14
medium, one in each direction. For > 1, partial re ections from the density discontinuity at the surface of the disk produce a series of Alfven wavefronts with successively decreasing amplitudes. These propagate into the external medium and cause an average exponential decline in the angular velocity of the disk. MP give analytic solutions valid for all values of . They examined two types of initial conditions: discontinuous initial conditions (DIC) in which the disk's rotational velocity 0 is uniform but drops discontinuously to zero at the disk surface, and continuous initial conditions (CIC) in which the disk is given a di erential rotational velocity which follows a cosine function and smoothly drops to zero. We examine one of each of these problems. To complete the speci cation of the DIC problem, we adopt 0 = 1 and a disk ten times denser than the external medium ( = 10). For the CIC problem we take = 10 again, but now the initial angular velocity in the disk is given by (10)
(z Z ) = 12 0[1 + cos(z=Z )]; with 0 = 1, and the angular velocity is assumed to vanish in the external medium (corresponding to = 1 in MP). For the DIC, the emergent waveform is a series of declining steps. The values of the rst few of these steps in angular velocity are provided in Table 3 (note that B = 0v for this model), along with their times of emergence from the disk. The analytic solution for the CIC problem is more involved. We refer the reader to MP for an in-depth description of the solution in each case, as well as for all of the analytic formulae needed to generate the two solutions. For the DIC problem, these are eqns. 27, 28, 33 in MP, while for the CIC problem the relevant equations are 39, and 42-58. Note that although a very large number of equations are needed to specify the analytic solution, particularly the CIC problem, we found only one typographical error. Eqn. [56i] in MP should read M4 = sin(= )sin( = ). For reference we provide FORTRAN subroutines to compute both the DIC and CIC solutions as an appendix to this paper. In Figures 7 and 8 we plot the numerical solutions generated by ZEUS-2D along with 15
the analytic solution generated by MP (plotted as a solid line). The numerical solution is computed using 300 grid points on the domain x 2 (0; 15). The DIC solution is shown at t = 13, immediately before the rst wavefront leaves the computational domain. The CIC solution is depicted at t = 15 just as the rst wavefront leaves the domain. When plotted as points overlaid on analytic curves the numerical results show good agreement with the analytic solutions. For the DIC, the wavefronts are resolved with about fourteen grid points; this is consistent with the performance of the second-order van Leer interpolation method for advection problems. We have found the discontinuities can be resolved with fewer grid points by using a higher-order interpolation method, although in this case small amplitude ringing occurs behind the wavefronts. For the CIC problem, close inspection reveals that the largest discrepancies between the numerical and analytic solutions occur at local extrema. This is demonstrated more clearly by plotting the errors themselves, as below. In Figures 9 and 10 we plot the absolute errors in the numerical solutions for the DIC and CIC problems using resolutions of both 75 and 600 grid points. For the DIC, the error is dominated by the nite resolution at discontinuities, as expected. The amplitude of the error is comparable for both resolutions, but the physical width of the error spike is smaller at the higher resolution; it spans the same number of grid points in both cases. For the CIC, the errors are more uniform, although the maximum error still occurs at the wavefronts. In this case, however, the error at higher resolution is uniformly lower by over two orders of magnitude than that at lower resolution. We also note that for the CIC, the errors do not increase at the outer boundary indicating that no re ection of the outgoing Alfven waves is occurring. As we have argued previously, measuring the convergence rate for numerical schemes in problems which contain discontinuities will yield a value of, at most, 1.0. Indeed, the convergence rate of the MOC-CT algorithm for the DIC problem is near one regardless of the order of interpolation used. We have, however, measured the convergence rate for the CIC problem using both the van Leer and PPA interpolation methods. The behavior of the 16
L1 error norm versus resolution for both schemes is given in Figure 11, along with a straight
line of slope -2 for comparison. The slightly larger error in the magnetic eld visible for both schemes results from the kink in the solution for the magnetic eld which occurs at the disk surface (see Figure 8). At a resolution of 600 zones the convergence rate is 2.1 for the van Leer scheme and 2.5 for the PPA scheme. The principle source of error in the PPA results appears to be small amplitude ringing which occurs immediately behind the discontinuities. 3.4
One-Dimensional Magnetic Stellar Wind Solutions
One important area of astrophysical application for MHD codes is to problems involving accretion and winds. Because of the interest in magnetic out ows there are many analytic studies available. We consider two such magnetized out ow problems in this paper. The rst of these is the steady 1D magnetized wind solution of Weber and Davis (1967), which assumes the poloidal magnetic eld geometry is purely radial. In general, for a wind that possesses angular momentum, cross- eld force balance will not admit a radial poloidal eld. However, suciently near the equatorial plane this can be a useful approximation. We follow the procedure outlined by Sakurai (1985) to obtain a speci c Weber-Davis wind model for use as a test problem. Sakurai (1985) has also described a method for obtaining twodimensional axisymmetric wind solutions, but these are dicult to obtain and we choose to con ne attention here to a wind with strictly radial poloidal eld. The assumption of steady one-dimensional ow allows a set of rst integrals to be obtained from the MHD equations that express conservation of speci c entropy, mass, magnetic
ux, ux-freezing, angular momentum, and energy. Following Sakurai these can be written as p = K ; (11) vr r2 = f;
(12)
Br r2 = 8;
(13)
17
(v 0 r)Br = vrB ; ! Br B 2 r v 0 4vr = rA ;
(14) (15)
1 [v2 + (v 0 r)2] + p 0 GM 0 1 2r2 = E: (16) 2 r
01 r 2 These contain six integration constants: K , f , 8, , rA, and E , not all of which are independent for a critical ow. Substituting equations [11]-[15] into equation [16] yields the Bernoulli equation H (r; ) = E that determines the density pro le (r) of the ow. A critical (wind) solution must simultaneously pass through both the slow and fast magnetosonic critical points and be smooth across the Alfven point (where the radial ow velocity equals the Alfven velocity). For such a wind, it is sucient to prescribe two constants, rA and A, the radius and density at the Alfven point, and two out of the following three dimensionless parameters 82 ; (17) = 4GMA r3 2=
A
01
KA rA
GM
2rA3 : != GM
;
(18)
(19) Here is the ratio of magnetic to gravitational potential energy, 2 is the ratio of thermal to gravitational potential energy and ! is the ratio of rotational to gravitational potential energy. We take 2 and ! as the speci able parameters. Then the four parameters (rA, A, 2, and !) and the two critical point conditions determine the remaining two constants, and E . Using these values in equations [17]-[19] and using the relation A = 4f 2=82, the constants that appear in [11]-[16], K , f , 8, and , can be found. By adopting the dimensionless radius, x = r=rA , and dimensionless density, y = =A , the Bernoulli equation can be written in dimensionless form [Sakurai 1985 eqn. (10)] ! 2 ( x 0 1 =x ) 1 ! 2 2
0 1 ~ (20) H~ = 4 2 + 2x y 0 1 y 0 x + 2 (y 0 1)2 0 x = E; ~ = @ H=@y ~ = 0, where E~ = ErA =GM . A wind solution is obtained by requiring that @ H=@x and H~ 0 E~ = 0 hold at the both the fast and slow critical points. Simultaneous solution 18
of these six equations (given 2 and !) determines the radii, xs and xf , and densities, ys and yf , of the slow and fast points, and and E~ . These nonlinear simultaneous equations can be solved using standard iterative techniques such as Newton-Raphson (see Press et al. 1986). Once all the parameters have been determined, the density pro le (r) can be obtained from the Bernoulli equation, also via an iterative procedure. Sucient accuracy is required in solving the critical point conditions to ensure that the pro le is smooth and continuous through the critical points. In order to create a speci c test problem, we adopt 2 = 0:5, ! = 0:25, and = 1:2 [one of the parameter sets used by Sakurai (1985)]. Solving the critical point conditions with these values gives a dimensionless energy E~ = 1:7384 and = 0:5756. The dimensionless radius and density (x; y) at the critical points are found to be (0:7770; 1:9405) for the slow magnetosonic point and (1:3024; 0:5139) for the fast. If we further specify that the dimensions of the problem be set by adopting rA = 500, A = 10013 and GM = 1, the remaining parameters can be determined. These are summarized in Table 4. An analytic solution such as the Weber-Davis wind can be exploited in several ways. The most obvious is to establish an in ow inner boundary condition and allow the code to evolve to the steady state, after which the numerical errors can be measured. An alternative procedure is to initialize the entire grid with the analytic solution. As this solution is evolved numerically, di erences will develop between the numerical and analytic solution. These di erences provide a direct measure of the truncation error. In fact, by placing the analytic solution on the grid initially, the numerical error can be analyzed in a single timestep. A further advantage of this test is that individual MHD equations can be separately analyzed in detail by keeping some terms xed at their analytic values and evolving others. For example one can x the velocity and study the errors that arise solely from density transport. Single timestep convergence to an analytic solution has proven to be a very powerful technique for uncovering bugs and algorithmic inadequacies. 19
Given the parameters for the problem in Table 4, the computational domain can be xed by taking an inner boundary at rin = 0:58rA , an outer boundary at 1:83rA and 100 grid zones equally spaced in radius. This places all three critical points well within the grid. We have run the computation to 2:84 2 104, long enough to establish a steady state. In Figure 12 we present the Weber-Davis test problem in the traditional manner: analytic curves with numerical values overlaid. Viewed in this way the code is producing \good" results and the errors are less than a percent in every variable. Figure 13 presents a more informative story by plotting the relative error for each variable as a function of radius. The typical error is about a tenth of a percent. Further, from these plots it is clear that the dominant error is due to the boundary conditions. When we increase the resolution and repeat this test problem we nd that the error everywhere goes as 1r, i.e., a rst-order rate of convergence. We next exploit the Weber-Davis problem to study further the sources of error. As a test, initialize the grid with the analytic solution and suspend the numerical computation of the velocity. We nd that the error in the density is reduced substantially, and the error decreases with resolution as 1r2 (using the velocity corrected transport algorithm for Eulerian transport; see Finn & Hawley 1989). This test demonstrates that the dominant error in the density comes from the computation of the velocity, as well as the fact that the numerical solution for any given variable (here density) can be no better than the weakest component. A more detailed examination of the toroidal eld also provides similar insight. Again, xing the velocity to have the analytic dependence we compute the evolution of B . The MOC procedure produces errors that converge at a rst-order rate, indicating that our present implementation is rst-order. This could be due solely to the boundary condition or it could be the order of the scheme. The question can be addressed through single-timestep convergence testing. In such a test most grid zones will not interact with the boundary and their error will be representative of the algorithm. Performing this test con rms the hypothesis that our MOC procedure as implemented is rst-order. In the Weber-Davis test 20
both the spherical geometry and the spatially-varying Alfven velocity introduce greater complexity into the equations. A possible area for further algorithmic re nement includes more accurate treatment of the spherical geometry terms ( , Monchmeyer & Muller 1989), and accounting for the subzone structure in the Alfven velocity in the characteristic calculation. e.g.
We perform an additional test by replacing the MOC scheme with a simple, formally second-order spatially-centered di erence procedure for the radial derivative of the angular velocity, @r . The Eulerian transport term in the B evolution equation is unchanged. If we perform a test with the velocity terms kept at their analytic values we nd errors that converge as 1r2. However, the use of spatially-centered di erencing for the Alfven terms is unacceptable in fully dynamic calculations because the highly dispersive nature of the resulting algorithm yields excessive high-wavenumber ringing in the presence of discontinuities. We note in passing that these high-wavenumber components are not immediately obvious in plots such as Figure 12 since the amplitude of the oscillations is small compared to B . They are more than obvious in an error plot such as Figure 13. 3.5
Multi-Dimensional Solar Coronal Transient Models
The nal problem in our test suite is a two-dimensional dynamical out ow that involves shocks. Multidimensional tests are especially important for verifying an MHD algorithm's performance in the presence of eld curvature, but nding suitable test problems is dicult. Nonetheless, time-dependent two- and three-dimensional models of solar coronal transients have been obtained by Low (1984, and references therein) and these serve as excellent test problems. Low's (1984) models were originally constructed to study solar coronal transients, the violent expulsions of coronal plasma that have been the subject of intense observational study in the last decade (e.g., Low 1990). The models are constructed by assuming that at t = t a hot magnetized bubble is ejected and forced into a spherically symmetric, non21
magnetic ambient medium. At the contact surface between the two media a strong shock develops and runs ahead of the contact surface, heating and compressing the non-magnetic ambient gas. Let RS be the radial position of the shock and RC be the position of the contact surface. Then at late times, the ow is characterized by three regions. For r < RC there is a spherically symmetric radial out ow and an axially symmetric magnetic eld. Stresses in this eld are balanced by the gravitational force and gradients in the poloidal pressure. The region outside of the magnetized uid, Rc < r < RS , consists of a spherically symmetric radial out ow composed of a shell of hot, compressed unmagnetized ambient gas. Beyond that at r > RS there is a spherically symmetric undisturbed, non-magnetic ambient medium. This problem tests the full radial dynamics and angular magnetic and pressure gradient force terms. Note that while we ran this test problem with a spherical-polar code, it is an even more demanding test problem in cylindrical coordinates. The spherically symmetric radial out ow in the region r < RS is generated by using the strong shock jump conditions and the properties of isentropic ow. From Low (1984), the radial out ow can be written as 8 > > p (r; t) = (4 3 )01 [GM=r]4 > < s (21) F or r < Rc > s (r; t) = 03 [GM=r]3 > > : vr (r; t) = r=t
where
8 > > p (r; t) = (7=6)804 d R2 exp[2GMR6 =3 9 ] > < s F or RC < r < RS > s (r; t) = 7803 d (R= )8 exp[2GMR6 =3 9 ] > > : vr (r; t) = r=t 8 > < s (r) = d (R =r)26=7 exp[(2GM=3R3 )(R =r)9=7 )] F or r > RS > : vr (r) = 0
8 = 8(t) = 1=2t; = (t) = R 81=6; 22
(22) (23)
RC RS d
= RC (t) = C 8; = RS (t) = R87=6; = 108mp exp[02GM=3R3 ]:
(24)
The axially symmetric magnetic medium participating in the out ow for r < RC is generated using the methods outlined by Low (1984). Thus, Br (r; ; t) B (r; ; t) B (r; ; t) p(r; ; t) (r; ; t)
where
= = = = =
2Ar02[P + ( )1=2J3=2( )]cos 0A(r8)01 [( )1=2J1=2( ) 0 ( )01=2J3=2( )]sin A(r8)01 [P + ( )1=2 J3=2( )]sin (4r4)01AP(2 0 2 2)A~(; ) + ps (r; t) (2GMr3 )01AP(4 0 2 2)A~(; ) + s (r; t) A~(; ) = A [P + ( )1=2 J3=2( )]sin2 ;
(25)
(26)
J1=2 and J3=2 are the modi ed Bessel functions of half-integer order, and ps and s are given
in eq. (21) above. Numerical values for all the parameters needed in the above analytic expressions are given in Table 5. The above expressions correct several minor typographical errors in the formulae given in Low (1984). To run the Low model as a test, we initialize the analytic solution given in eqs. (2126) with t = t0 = 87400 on a grid in spherical polar coordinates using 200 radial by 120 angular grid points on the domain r 2 (1 2 1011; 5 2 1011) and 2 (0; ). Note that although Low's models possess equatorial symmetry, we evolve the model over the full angular domain to ensure the algorithm can maintain this symmetry. It is also useful to run the problem on the angular domain 2 (0; =2) to ensure the symmetry boundary condition at the equatorial plane is implemented properly. Boundary conditions at the lower radial boundary are speci ed using the time dependent analytic solution, while the outer radial boundary is an out ow condition. The polar angular boundaries are both axes of symmetry. The initial 23
conditions are evolved to a time of t = 2:62 2 104 s, which coincides with the time of the nal plot given by Low (1984; gs. 5-7). In addition to comparing our results with these plots, we can compute the numerical error using the analytic solution. Figure 14 shows the density, poloidal eld lines and toroidal eld component at three di erent times during the evolution. Qualitatively, the numerical solution agrees well with Low (1984); the ejection of a magnetic bubble can be clearly followed. In Figure 15 we present the relative errors in the density, radial velocity, and radial and toroidal components of the magnetic eld as compared to the analytic solution at the termination of the calculation. We note a disadvantage in using the relative error as previously de ned since it diverges at the zeroes of the analytic functions. We nd the errors are dominated by two features. A large, well resolved wave of error propagates inward; it is generated by the inability of the numerical solution to represent properly the discontinuities present in the solution at t = t . In addition, there is a small, unresolved ring of error that propagates inward located roughly at r = 3 2 1011 which was generated by the partial re ection of the outgoing shock at the outer boundary. There is also a large error in each variable at the outer boundary. Both of these latter features are an indication of the limitations in the ability of our out ow boundary condition to transmit compressive waves. Improvements might be needed in some applications. Figure 16 is the L1 error norm for each variable as a function of resolution along the radial coordinate. We have scaled the errors by = 1:67 2 10016 gm cm03, p = 3:49 2 1003 dynes, B = 0:3 G, and v = 1:14 2 107 cm s01 so that all errors are roughly of the same order of magnitude. To provide a qualitative measure of the convergence rate, we also plot a straight line with a slope of -1. We note that all variables are converging at a rate of 1 for this problem.
24
4
Conclusions
In this paper we provide a comprehensive test suite for the calibration of MHD algorithms. This test suite should be useful for MHD code development as it provides a quantitative means of comparing di erent MHD algorithms. The tests include simple linear advection, an MHD shock tube, Alfven wave propagation, a magnetic wind, and a twodimensional coronal transient out ow. The linear advection test is a standard component of a test suite. It provides a measure of the numerical di usion associated with the advection algorithm. As is well known, there are large di erences between rst, second, and third order interpolation methods; only the higher-order methods are suitable for simulations. What makes eld advection especially important for MHD is the fact that gradients of B enter into the Lorentz force computation. It is therefore important to consider the current density pro les that are generated by the scheme. A constraint on the scheme should be that it not introduce anomalous current reversals near discontinuities, as occasionally occurs with the use of numerical \ ux correctors" or monotonicity constraints. The MHD Riemann problem of Brio and Wu (1988) provides an analogue to the Sod shock tube problem that is frequently used in calibrating hydrodynamical algorithms. The problem rigorously tests the propagation of strong compressive MHD waves. While the lack of an analytic solution poses some limit to the value of this test, we have shown that plots of relative error and measurements of the self-convergence rate provide a quantitative means of comparing and contrasting code results. The magnetic braking of an aligned rotator test problem measures the ability of an MHD scheme to propagate Alfven waves. This test focuses exclusively on Alfven waves and is likely to be a crucial gauge of the numerical method since the dynamics of many astrophysical systems is determined by the propagation of these waves ( the magnetic braking of an accretion disk). We found that some schemes calculate compressive MHD e.g.
25
waves accurately and provide good results for the MHD Riemann problem, but fail the Alfven wave test completely, pointing thus to the need for comprehensive testing. The magnetic stellar wind solution of Weber and Davis (1967) provides a test problem of considerable astrophysical relevance. This problem tests pressure, magnetic, centrifugal and gravitational acceleration terms, advection in spherical coordinates with nonconstant velocities, and toroidal Alfven waves. By various choices of grid size, problem parameters, and critical point locations one can tailor the test to emphasize speci c terms in the MHD equations. For example, increasing 2 increases the importance of thermal acceleration, while the importance of the magnetic eld can be emphasized by raising . Convergence testing, both to an extended time and with single timesteps, can be used to guarantee errorfree implementation of numerical algorithms. Measurements of error and convergence rates provide a quantitative means of algorithm comparison. The solar coronal transient solution of Low (1984) provides the only completely twodimensional problem in the test suite. By balancing eld tension forces against gravity and pressure gradients it tests the multi-dimensional eld curvature e ects that occur in MHD
ows. The existence of an analytic solution provides the same opportunities as the solar wind problem for code analysis and calibration. Taken together, the problems described in this paper provide a rather comprehensive test of MHD algorithms. Of course, these problems do not represent a complete set and continued testing of a numerical algorithm in concert with its use in applications is advised. For example, the self-convergence method outlined in x3.2 can be used during an application to determine that the changes in the numerical solutions are indeed decreasing with resolution, and the resolution is sucient to ensure that the physical e ects of interest are accurately determined. In this sense testing remains an integral and ongoing part of the process of numerical modeling. It is with great pleasure that we acknowledge discussions with B.C. Low, Dimitri Miha26
las and Telemachos Mouschovias. J. S. would like to thank Dimitri Mihalas and the National Center for Supercomputing Applications (NCSA) for nancial support. J. F. H. acknowledges support under NASA grant NAGW-1510 and NSF grant PHY90-18251. C. R. E. acknowledges support from NSF grant PHY90-01645, Presidential Young Investigator award PHY90-57865, and the Alfred P. Sloan Foundation. Simulations were performed on the Cray systems at the NCSA.
27
APPENDIX
As a convenience to other researchers, in this appendix we give FORTRAN listings of subroutines which will compute the solution to the magnetic braking of an aligned rotator as given by Mouschovias and Paleologou (1980) for both discontinuous and continuous initial conditions. c----------------------------------------------------------------------c SUBROUTINE DIC c subroutine dic(t,rho,o0,r,is,ie,x1b,omega,bphi) c c Computes the solution to the aligned magnetic braking problem c with discontinuous initial conditions (DIC) as given by Mouschovias c and Paleologou (1980). Equation numbers and variable names refer to c that paper. c--Input Variables: c t = dimensionless time (tau) c rho = dimensionless density c o0 = (omega)o - initial angular velocity of disk c r = dimensionless radius (zeta) c is,ie = starting and ending array indices at which soln is computed c x1b(i) = axial (z) positions of solution points c--Output Variables: c omega(i) = analytic solution for angular vel. at every grid point c bphi(i) = analytic solution for phi mag. fld. at every grid point c implicit undefined(a-z) integer in parameter (in=100) c integer is,ie real t,rho,o0,r,x1b(in),omega(in),bphi(in) c integer i,n,m real sr,z,l1,l2,l3 c======================================================================= sr = sqrt(rho) c c do the case rho > 1.0 c if (rho .gt. 1.0) then do 100 i=is,ie z = x1b(i) if (z .le. 1.0) then if (t .lt. sr*(1.0-z)) then omega(i) = o0 bphi (i) = 0.0 else
28
10
20
30
40
100
n = 0 continue l1 = sr*(2.0*float(n) + 1.0 - z) l2 = sr*(2.0*float(n) + 1.0 + z) l3 = sr*(2.0*float(n) + 3.0 - z) if (t .ge. l1 .and. t .lt. l2) then omega(i) = o0*sr/(sr+1.0)*((sr-1.0)/(sr+1.0))**n bphi (i) = -r*o0*sr/(sr+1.0)*((sr-1.0)/(sr+1.0))**n goto 20 endif if (t .ge. l2 .and. t .lt. l3) then omega(i) = o0*((sr-1.0)/(sr+1.0))**(n+1) bphi (i) = 0.0 goto 20 endif n = n+1 goto 10 continue endif else if (t .lt. z-1.0) then omega(i) = 0.0 else m = 0 continue l1 = 2.0*float(m )*sr + z - 1.0 l2 = 2.0*float(m+1)*sr + z - 1.0 if (t .ge. l1 .and. t .lt. l2) then omega(i) = o0*sr/(sr+1.0)*((sr-1.0)/(sr+1.0))**m goto 40 endif m = m+1 goto 30 continue endif bphi(i) = -r*omega(i) endif continue endif
c c do the case rho = 1.0 c if (rho .eq. 1.0) then do 200 i=is,ie z=x1b(i) if (z .le. 1.0) then if (t .lt. 1.0-z) then omega(i) = o0 bphi (i) = 0.0 endif if (t .ge. 1.0-z .and. t .lt. 1.0+z) then omega(i) = o0/2.0 bphi (i) = -r*o0/2.0 endif
29
200 c
if (t .ge. omega(i) bphi (i) endif else if (t .lt. omega(i) bphi (i) endif if (t .ge. omega(i) bphi (i) endif if (t .ge. omega(i) bphi (i) endif endif continue endif
1.0+z) then = 0.0 = 0.0 z-1.0) then = 0.0 = 0.0 z-1.0 .and. t .lt. z+1.0) then = o0/2.0 = -r*o0/2.0 z+1.0) then = 0.0 = 0.0
return end c----------------------------------------------------------------------c SUBROUTINE CIC c subroutine cic(t,rho,o0,r,is,is,x1b,beta,omega,bphi) c c Computes the solution to the aligned magnetic braking problem c with continuous initial conditions (CIC) as given by Mouschovias c and Paleologou (1980). Equation numbers and variable names refer to c that paper. c--Input Variables: c t = dimensionless time (tau) c rho = dimensionless density c o0 = (omega)o - initial angular velocity of disk c r = dimensionless radius (zeta) c is,ie = starting and ending array indices at which soln is computed c x1b(i) = axial (z) positions of solution points c beta = smoothing length of initial angular velocity c--Output Variables: c omega(i) = analytic solution for angular vel. at every grid point c bphi(i) = analytic solution for phi mag. fld. at every grid point c implicit undefined(a-z) integer in real pi parameter (in=100,pi=3.141592654) c integer is,ie real t,rho,o0,r,x1b(in),beta,omega(in),bphi(in) c integer i real sr,z,i1,i2,i3,i4,i5,j1,j2,j3,j4,k1,k2,k3,l5
30
external geti,getj,getk c======================================================================= sr = sqrt(rho) c c do the case rho > 1.0 c if (rho .gt. 1.0) then do 100 i=is,ie z = x1b(i) if (z .le. 1.0) then call geti(t,rho,z,beta,i1,i2,i3,i4,i5) l5 = -sin(pi*z/beta)*sin(pi*t/(sr*beta)) omega(i) = 0.5*o0*(i1+i2+i3+i4+i5) bphi (i) = r*sr*0.5*o0*(i1-i2+i3-i4+l5) endif if (z .gt. 1.0 .and. z .le. beta) then call getj(t,rho,z,beta,j1,j2,j3,j4) l5 = -sin(pi*z/beta)*sin(pi*t/beta) omega(i) = 0.5*o0*(j1+j2+j3+j4) bphi (i) = r*0.5*o0*(j1-j2-j3+l5) endif if (z .gt. beta) then call getk(t,rho,z,beta,k1,k2,k3) omega(i) = 0.5*o0*(k1+k2+k3) bphi (i) = -r*omega(i) endif 100 continue endif c c do the case rho = 1.0 c if (rho .eq. 1.0) then do 200 i=is,ie z=x1b(i) if (z .le. beta) then if (t .lt. beta-z) then omega(i) = 0.5*o0*(1.0+cos(pi*z/beta)*cos(pi*t/beta)) bphi (i) = -r*0.5*o0*sin(pi*z/beta)*sin(pi*t/beta) endif if (t .ge. beta-z .and. t .lt. beta+z) then omega(i) = 0.5*o0*(cos(pi*(t-z)/(2.0*beta)))**2 bphi (i) = -r*0.5*o0*(cos(pi*(t-z)/(2.0*beta)))**2 endif if (t .ge. beta+z) then omega(i) = 0.0 bphi (i) = 0.0 endif else if (t .lt. z-beta) then omega(i) = 0.0 bphi (i) = 0.0 endif if (t .ge. z-beta .and. t .lt. z+beta) then omega(i) = 0.5*o0*(cos(pi*(t-z)/(2.0*beta)))**2
31
200 c
bphi (i) endif if (t .ge. omega(i) bphi (i) endif endif continue endif
= -r*0.5*o0*(cos(pi*(t-z)/(2.0*beta)))**2 z+beta) then = 0.0 = 0.0
return end c----------------------------------------------------------------------c SUBROUTINE GETI c subroutine geti(t,rho,z,beta,i1,i2,i3,i4,i5) c c Computes the Ik integrals defined by eqns 42, and 44-47 in MP c implicit undefined(a-z) real pi parameter(pi=3.141592654) real t,rho,z,beta,i1,i2,i3,i4,i5 c integer ii1,ii2,ii3,ii4 real sr,bd,br,x1,x2,x3,x4,l1,l2,l3 c integer ii real dd,pp,rr,xx real eqn42b,eqn45b eqn42b(dd,pp,rr,xx,ii) = . cos(pp*pi/beta)/dd*(cos(pi*(xx+1.0)/beta) . - rr**ii*cos(pi*(xx+1.0-2.0*pp*float(ii))/beta)) . + pp*sin(pp*pi/beta)/dd*(sin(pi*(xx+1.0)/beta) . - rr**ii*sin(pi*(xx+1.0-2.0*pp*float(ii))/beta)) . - 0.5*cos(pi*xx/(pp*beta)) . + 0.5*rr**ii*cos(pi*(xx-2.0*pp*float(ii))/(pp*beta)) eqn45b(dd,pp,rr,xx,ii) = . cos(pp*pi/beta)/dd*(cos(pi*xx/beta) . - rr**ii*cos(pi*(xx-2.0*pp*float(ii))/beta)) . + pp*sin(pp*pi/beta)/dd*(sin(pi*xx/beta) . - rr**ii*sin(pi*(xx-2.0*pp*float(ii))/beta)) . + 0.5*(rr**ii - 1.0) c sr = sqrt(rho) bd = 2.0*((cos(sr*pi/beta))**2 + rho*(sin(sr*pi/beta))**2) br = (sr-1.0)/(sr+1.0) c compute i1 integral if (t .lt. sr*(1.0-z)) then i1 = 0.0 else x1 = t + sr*z ii1 = 1 10 continue
32
20 c
30
40 c
50
60 c
70
l1 = t-sr*(1.0-z) l2 = 2.0*sr*float(ii1) l3 = t+sr*(1.0+z) if (l2 .gt. l1 .and. l2 .le. l3) goto ii1 = ii1 + 1 goto 10 continue i1 = eqn42b(bd,sr,br,x1,ii1) endif compute i2 integral if (t .lt. sr*(1.0+z)) then i2 = 0.0 else x2 = t - sr*z ii2 = 1 continue l1 = t-sr*(1.0+z) l2 = 2.0*sr*float(ii2) l3 = t+sr*(1.0-z) if (l2 .gt. l1 .and. l2 .le. l3) goto ii2 = ii2 + 1 goto 30 continue i2 = eqn42b(bd,sr,br,x2,ii2) endif compute i3 integral if (t .lt. sr*(1.0-z)+(beta-1.0)) then i3 = 0.0 else x3 = t + sr*z - beta+1.0 ii3 = 1 continue l1 = t - sr*(1.0-z) - (beta-1.0) l2 = 2.0*sr*float(ii3) l3 = t + sr*(1.0+z) - (beta-1.0) if (l2 .gt. l1 .and. l2 .le. l3) goto ii3 = ii3+1 goto 50 continue i3 = eqn45b(bd,sr,br,x3,ii3) endif compute i4integral if (t .lt. sr*(1.0+z)+(beta-1.0)) then i4 = 0.0 else x4 = t - sr*z - beta+1.0 ii4 = 1 continue l1 = t - sr*(1.0+z) - (beta-1.0) l2 = 2.0*sr*float(ii4) l3 = t + sr*(1.0-z) - (beta-1.0) if (l2 .gt. l1 .and. l2 .le. l3) goto ii4 = ii4+1 goto 70
33
20
40
60
80
80
continue i4 = eqn45b(bd,sr,br,x4,ii4) endif c compute i5 integral i5 = 1.0 + cos(pi*z/beta)*cos(pi*t/(sr*beta)) c return end c----------------------------------------------------------------------c SUBROUTINE GETJ c subroutine getj(t,rho,z,beta,j1,j2,j3,j4) c c Computes the Jk integrals defined by eqns 48,49,51, and 52 in MP c implicit undefined(a-z) real pi parameter(pi=3.141592654) real t,rho,z,beta,j1,j2,j3,j4 c integer jj real sr,dd,rr,yy,ff,xx,l1,l2,l3 c sr = sqrt(rho) dd = 2.0*((cos(sr*pi/beta))**2 + rho*(sin(sr*pi/beta))**2) rr = (sr-1.0)/(sr+1.0) ff = 2.0*sr/(sr+1.0) c compute j1 integral if (t .lt. beta-z) then j1 = 0.0 else j1 = -1.0*(cos(pi*(t+z)/(2.0*beta)))**2 endif c compute j2 integral if (t .lt. z-1.0) then j2 = 0.0 else yy = t-z+1.0 jj = 1 10 continue l1 = t-z+1.0 l2 = 2.0*sr*float(jj) l3 = t-z+1.0+2.0*sr if (l2 .gt. l1 .and. l2 .le. l3) goto 20 jj = jj+1 goto 10 20 continue xx = pi/beta j2 = -2.0/dd*(sin(xx)*cos(sr*xx) - sr*cos(xx)*sin(sr*xx)) . *(cos(sr*xx)*sin(yy*xx) - sr*sin(sr*xx)*cos(yy*xx)) . - rr**(jj-1)*ff*( . (sr*sin(xx)*sin(sr*xx) . + cos(xx)*cos(sr*xx))*cos(xx*(yy-2.0*sr*float(jj)+sr))/dd . + (sr*cos(xx)*sin(sr*xx)
34
c
30
40
c c
. - sin(xx)*cos(sr*xx))*sin(xx*(yy-2.0*sr*float(jj)+sr))/dd . - 0.5*cos(xx*(yy/sr-2.0*float(jj)+1.0))) endif compute j3 integral if (t .lt. z+beta-2.0) then j3 = 0.0 else yy = t-z-beta+2.0 jj = 1 continue l1 = t-z-beta+2.0 l2 = 2.0*sr*float(jj) l3 = t-z-beta+2.0+2.0*sr if (l2 .gt. l1 .and. l2 .le. l3) goto 40 jj = jj+1 goto 30 continue xx = pi/beta j3 = (((cos(sr*xx))**2-rho*(sin(sr*xx))**2)*cos(xx*yy) . + 2.0*sr*sin(sr*xx)*cos(sr*xx)*sin(xx*yy))/dd . - 0.5 + rr**(jj-1)*ff* . (0.5 - sr/dd*sin(sr*xx)*sin(xx*(yy-2.0*sr*float(jj)+sr)) . - cos(sr*xx)/dd*cos(xx*(yy-2.0*sr*float(jj)+sr))) endif compute j4 integral j4 = 1.0 + cos(pi*z/beta)*cos(pi*t/beta)
return end c----------------------------------------------------------------------c SUBROUTINE GETK c subroutine getk(t,rho,z,beta,k1,k2,k3) c c Computes the Kk integrals defined by eqns 49,51, and 53 in MP c implicit undefined(a-z) real pi parameter(pi=3.141592654) real t,rho,z,beta,k1,k2,k3 c integer jj real sr,dd,rr,yy,ff,xx,l1,l2,l3 c sr = sqrt(rho) dd = 2.0*((cos(sr*pi/beta))**2 + rho*(sin(sr*pi/beta))**2) rr = (sr-1.0)/(sr+1.0) ff = 2.0*sr/(sr+1.0) c compute k1 integral (identical to j2 integral) if (t .lt. z-1.0) then k1 = 0.0 else yy = t-z+1.0 jj = 1
35
10
20
c
30
40
c
c
continue l1 = t-z+1.0 l2 = 2.0*sr*float(jj) l3 = t-z+1.0+2.0*sr if (l2 .gt. l1 .and. l2 .le. l3) goto 20 jj = jj+1 goto 10 continue xx = pi/beta k1 = -2.0/dd*(sin(xx)*cos(sr*xx) - sr*cos(xx)*sin(sr*xx)) . *(cos(sr*xx)*sin(yy*xx) - sr*sin(sr*xx)*cos(yy*xx)) . - rr**(jj-1)*ff*( . (sr*sin(xx)*sin(sr*xx) . + cos(xx)*cos(sr*xx))*cos(xx*(yy-2.0*sr*float(jj)+sr))/dd . + (sr*cos(xx)*sin(sr*xx) . - sin(xx)*cos(sr*xx))*sin(xx*(yy-2.0*sr*float(jj)+sr))/dd . - 0.5*cos(xx*(yy/sr-2.0*float(jj)+1.0))) endif compute k2 integral (identical to j3 integral) if (t .lt. z+beta-2.0) then k2 = 0.0 else yy = t-z-beta+2.0 jj = 1 continue l1 = t-z-beta+2.0 l2 = 2.0*sr*float(jj) l3 = t-z-beta+2.0+2.0*sr if (l2 .gt. l1 .and. l2 .le. l3) goto 40 jj = jj+1 goto 30 continue xx = pi/beta k2 = (((cos(sr*xx))**2-rho*(sin(sr*xx))**2)*cos(xx*yy) . + 2.0*sr*sin(sr*xx)*cos(sr*xx)*sin(xx*yy))/dd . - 0.5 + rr**(jj-1)*ff* . (0.5 - sr/dd*sin(sr*xx)*sin(xx*(yy-2.0*sr*float(jj)+sr)) . - cos(sr*xx)/dd*cos(xx*(yy-2.0*sr*float(jj)+sr))) endif compute k3 integral if (t .lt. z-beta) then k3 = 0.0 else k3 = (sin(pi*(t-z+beta)/(2.0*beta)))**2 endif return end
36
TABLE 1
Variable Left 1.0 P 1.0 vx 0 vy 0 By 1.0
MHD Shock Tube Values FR SC CDl CDr 0.664 1 0.817 0.701 0.240 (0.676) (0.697) 0.443 0.687 0.509 0.509 (0.457) (0.516) 0.662 0.484 0.597 0.597 (0.637) (0.599) -0.248 -1.16 -1.58 -1.58 (-0.233) (-1.58) 0.567 -0.179 -0.536 -0.536 (0.585) (-0.534)
FR Right 0.116 .125 0.089 .1 -0.255 0 -0.179 0 -0.896 -1.0
1 values in brackets are those given by Brio and Wu (1988).
TABLE 2
MHD Shock Tube Density Errors No. Zones Integrated Error Convergence Rate 64 4.1137E-03 .... 128 1.9829E-03 1.04 256 9.7396E-04 1.02 512 4.8265E-04 1.01 1024 2.4014E-04 1.01 TABLE 3
Magnetic Braking DIC Solution v Emergence Time 0.7598 0.00 0.3947 6.33 0.2050 12.65 0.1065 18.97 0.0553 25.30 0.0288 31.62 37
TABLE 4
Weber-Davis Wind Parameters Parameter Value Parameter 2 0.5 GM ! 0.25
1.2 K 0.5756 8 E~ 1.7384 f rA 500 E A 10013 vA
Value 1.0 4:472 2 1005 0.3318 9:509 2 1003 8:482 2 10010 3:477 2 1003 3:393 2 1002
TABLE 5
Values for parameters in Low model (in cgs units). 5:54 2 1011 cm01 A 1:5 2 1021 Gauss cm2 P 1:01327 C 1:104 2 1011 cm 5:24 2 1008 sec02 2:42 2 1020 cm3sec02gm01=3 M 2:0 2 1033 gm R 1:0 2 1011 cm t 8:74 2 103 sec
38
REFERENCES
Beck, R, Kronberg, P., & Wielebinski, R. (eds.) 1990, Galactic and Extragalactic Magnetic Fields, (Dordrecht: Riedel) Belvedere, G. (ed) 1989, Accretion Disks and Magnetic Fields in Astrophysics, (Dordrecht: Kluwer) Blandford, R. B. 1989, in Theory of Accretion Disks, ed. F. Meyer, W. J. Duschl, J. Frank, and E. Meyer-Hofmeister, (Dordrecht: Kluwer), p. 35 Brio, M., & Wu, C.C. 1988, J Comp Phys, 75, 400 (BW) Clarke, D.A., Norman, M.L., & Burns, J.O. 1986, ApJ, 311, L63 Clarke, D.A., Norman, M.L., & Burns, J.O. 1989, ApJ, 342, 700 Colella, P. & Woodward, P., 1984, J Comp Phys, 54, 174 Courant, R., & Friedrichs, K.O. 1948, Supersonic Flow and Shock Waves, (New York: Springer-Verlag) Evans, C.R., & Hawley, J.F. 1988, ApJ, 33, 659 Finn, L. S., & Hawley, J. F. 1989, unpublished report Hawley, J. F., Smarr, L. L., and Wilson, J. R. 1984, ApJ Supp, 55, 211 Low, B.C. 1984, ApJ, 281, 392 Low, B.C. 1990, Ann. Rev. Astron. Astrophys., 28, 491 Monchmeyer R. & Muller E. 1989, A&A, 217, 351 Mouschovias, T.Ch. 1987, in Physical Processes in Interstellar Clouds, eds. G. Mor ll & M. Scholer, (Dordrecht:Riedel), 453 Mouschovias, T.Ch., & Paleologou, E.V. 1980, ApJ, 237, 877 (MP) Norman, M.L., Stone, J.M., Evans, C.R., & Hawley, J.F. 1991, in preparation 39
Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T. 1986, Numerical Recipes, (New York: Cambridge) Sakurai, T. 1985, A&A, 152, 121 Shu, F.H., Adams, F.C. & Lizano, S. 1987, Ann. Rev. Astron. Astrophys., 25, 23 Sod, G.A. 1978, J Comp Phys, 27, 1 Stone, J.M., & Norman, M.L. 1991a, ApJ Supp, submitted Stone, J.M., & Norman, M.L. 1991b, ApJ Supp, submitted Stone, J.M., Mihalas, D., & Norman, M.L. 1991, ApJ Supp, submitted van Leer, B. 1977, J Comp Phys, 23, 276 Wardle, M., & Konigl, A. 1990, ApJ, 362, 120 Weber E. J., & Davis L. 1967, ApJ, 148, 217 Zachary, A., & Colella, P. 1990, J Comp Phys, in press
40
Figure Captions
Fig. 1. Results for the advection of a square pulse of transverse magnetic eld originally 50 zones wide a distance of ve times its width using the (a) donor cell, (b) van Leer, and (c) PPA interpolation methods. The analytic locations of the eld discontinuity are 0.4 and 0.6. Fig. 2. Current density obtained from the square eld pulse computations of Figure 1. Analytically, the current should be a delta function at x = 0:4. The sharpest current density distribution is obtained from the PPA computation. A slightly broader distribution is produced by the van Leer scheme. The donor cell advection scheme produces a current spread out over a distance roughly equal to the original pulse. Fig. 3. The relative error for the advection test pictured in Figure 1 as a function of spatial location x. The PPA scheme produces the lowest error, the van Leer scheme slightly larger errors, and the donor cell scheme produces substantial error over a large spatial domain. Fig. 4. Magnetic Shock Tube Results. The (a) density, (b) pressure, (c) x-component of velocity, (d) y-component of velocity, and (e) y-component of magnetic eld are shown. The density plot shows the location of the (1) Fast Rarefaction wave, (2) Slow Compound wave, (3) Contact Discontinuity, (4) Slow Shock, and (5) Fast Rarefaction wave. The plots show a calculation with 800 grid zones, a Courant number of 0.5, a grid spacing of 1x = 0:125, and a time of 10. Fig. 5. The L1 error norm vs. number of grid zones for density in the MHD Riemann problem as a function of log of the number of grid zones (resolution). The errors are computed from the di erence between one numerical solution and a numerical solution computed at twice the resolution. Fig. 6. Di erence between values of density in the MHD Riemann problem computed using resolutions of 200 (solid line) and 800 zones (dashed line) as compared to a solution computed using 1600 zones. Discontinuities produce sharp spike-like error peaks; errors at rarefaction waves are smooth. 41
Fig. 7. Results for the magnetic braking of an aligned rotator starting from DIC and using the van Leer interpolation method in the MOC-CT scheme. Both (a) the rotational velocity and (b) the toroidal component of the magnetic eld are shown. The density ratio of the disk (located at x 1) to the ambient medium is ten, and 300 grid points are used. The analytic solution given by MP is plotted as a solid line. Fig. 8. Results for the magnetic braking of an aligned rotator starting from CIC and using the van Leer interpolation method in the MOC-CT scheme. Both (a) the rotational velocity and (b) toroidal component of the magnetic eld are shown. The density ratio of the disk (located at x 1) to the ambient medium is ten, and 300 grid points are used. The analytic solution given by MP is plotted as a solid line. Fig. 9. Di erences between the numerical and analytic solutions for the magnetic braking of an aligned rotator starting from DIC using 75 (top line) and 600 (bottom line) grid points for (a) the angular velocity, and (b) the toroidal magnetic eld. Fig. 10. Di erences between the numerical and analytic solutions for the magnetic braking of an aligned rotator starting from CIC using 75 (top line) and 600 (bottom line) grid points for (a) the angular velocity, and (b) the toroidal magnetic eld. Fig. 11. The L1 error norms for toroidal velocity and magnetic eld as a function of number of zones (resolution) for the magnetic braking of an aligned rotator starting from CIC. The lines labeled VL correspond to MOC-CT using the van Leer interpolation scheme, while the lines labeled PPA use MOC-CT with PPA. A reference line with slope of -2 (corresponding to a convergence rate of two) is included as a dashed line. Fig. 12. Plots of (a) density, (b) radial velocity, (c) toroidal magnetic eld, and (d) angular velocity as a function of radius in the Weber-Davis wind problem. The solid lines are the analytic values while the circles represent the numerical results from a 100 grid zone computation. The Alfven critical point is located at r = 500. Fig. 13. Plots of relative error in (a) density, (b) radial velocity, (c) toroidal magnetic eld, 42
and (d) angular velocity as a function of radius for the Weber-Davis wind problem. Fig. 14. Numerical solution of the solar coronal transient models. Plots (a) through (c) show contours of the density (solid lines) and representative poloidal eld lines projected into the poloidal plane (dashed lines) at times (a) t = 8:74 2 103 s, (b) t = 1:75 2 104 s, and (c) t = 2:62 2 104 s. Plots (d) through (f) show contours of the toroidal component of the magnetic eld at the same times. Fig. 15. Relative errors in the (a) density, (b) radial velocity, (c) radial magnetic eld, and (d) toroidal magnetic eld in the numerical solution to the coronal transient model computed using 200 radial and 120 angular zones at a time t = 2:62 2 104 s. The maximum errors in each case are 0.98%, 0.98%, 2.6%, and 2.1%. Sixteen equally spaced contours are shown between the maximum and zero. Fig. 16. The L1 error norm as a function of number of zones (resolution) for each of the dependent variables in the solar coronal transient model. The curves are labeled e for internal energy, v1 for radial velocity, d for density, B3 for toroidal eld, B2 for theta eld and B1 for radial eld. A reference straight line with a slope of -1 (corresponding to a convergence rate of one) is included as a dashed line.
43
Authors' Postal Addresses
Charles R. Evans University of North Carolina Dept. of Physics and Astronomy Phillips Hall CB #3255 Chapel Hill, NC, 27599 John F. Hawley University of Virginia Dept. of Astronomy P.O. Box 3818, University Station Charlottesville, VA, 22903 Michael L. Norman and James M. Stone National Center for Supercomputing Applications 5600 Beckman Institute, Drawer 25 405 N. Mathews Ave. Urbana, IL, 61801
44
