DISCO: A 3D Moving-Mesh Magnetohydrodynamics Code Designed for the Study of Astrophysical Disks
aa r X i v : . [ phy s i c s . c o m p - ph ] J un Draft version June 27, 2016
Preprint typeset using L A TEX style emulateapj v. 04/17/13
DISCO: A 3D MOVING-MESH MAGNETOHYDRODYNAMICS CODE DESIGNED FOR THE STUDY OFASTROPHYSICAL DISKS
Paul C. Duffell
Astronomy Department and Theoretical Astrophysics Center, University of California, Berkeley
Draft version June 27, 2016
ABSTRACTThis work presents the publicly available moving-mesh magnetohydrodynamics code DISCO. DISCOis efficient and accurate at evolving orbital fluid motion in two and three dimensions, especially at highMach number. DISCO employs a moving-mesh approach utilizing a dynamic cylindrical mesh that canshear azimuthally to follow the orbital motion of the gas. The moving mesh removes diffusive advectionerrors and allows for longer timesteps than a static grid. Magnetohydrodynamics is implemented inDISCO using an HLLD Riemann solver and a novel constrained transport scheme which is compatiblewith the mesh motion. DISCO is tested against a wide variety of problems, which are designed to testits stability, accuracy and scalability. In addition, several magnetohydrodynamics tests are performedwhich demonstrate the accuracy and stability of the new constrained transport approach, includingtwo tests of the magneto-rotational instability (MRI); one testing the linear growth rate and the otherfollowing the instability into the fully turbulent regime.
Subject headings: hydrodynamics — accretion disks — planetary systems: protoplanetary disks —X-rays: binaries — black hole physics — methods: numerical INTRODUCTION
The study of gaseous disks is of fundamental impor-tance to astrophysics. Disks are ubiquitous; essentiallyany gaseous orbital system which efficiently loses energywhile conserving angular momentum will form a disk.The formation of planets takes place in a protoplanetarydisk, which can influence the first few million years of theplanets’ existence (e.g. Kley & Nelson 2012; Testi et al.2014; Alexander et al. 2014). These especially include“transition disks” whose cavities may be signposts ofplanet formation (e.g. Espaillat et al. 2014). Accretiondisks around black holes (in particular, X-Ray binaries)are the most robust means of a black hole’s detectionvia electromagnetic waves (e.g. Remillard & McClintock2006; Abramowicz & Fragile 2013). Some stars are alsothought to be surrounded by an accreting disk (e.g.Rivinius et al. 2013). Accretion disks are also an im-portant feature of cataclysmic variables (e.g. Robinson1976; Papaloizou & Pringle 1978), and they are thoughtto be the power source behind active galactic nuclei (e.g.Osterbrock 1993; Sulentic et al. 2000). Current effortsto observe the horizon of Sgr A* depend on an emit-ting disk of gas surrounding this supermassive black hole(e.g. Psaltis & Johannsen 2011). Circumbinary disks areamong the most promising possibilities for electromag-netic counterparts of gravitational wave emission frommerging black holes (e.g. D’Orazio et al. 2013; Rafikov2013; Farris et al. 2014; Abbott et al. 2016; Stone et al.2016). Much of the same circumbinary physics applies toa newly-born binary star system, surrounded by a com-mon protostellar disk (e.g. Monin et al. 2007). Galaxiescan be considered another important type of disk (e.g.Rix & Bovy 2013; Graham 2013). Even Saturn’s ringsconstitute a disk, though it is not composed of gas, butof icy solids (e.g. Goldreich & Tremaine 1982).Each particular instance of a disk in nature possesses duff[email protected] its own specific physical ingredients. Many disks are ion-ized and therefore magnetic fields are important to theirevolution. Protoplanetary disks are dusty and also sub-ject to the gravitational influence of the planets beingformed in the disk. Black hole accretion disks can con-stitute very extreme environments where radiation hy-drodynamics, weak sector couplings, and general relativ-ity all come into play. In galaxies, self-gravity is veryimportant, as opposed to many other systems which canbe approximated as orbiting a single point mass at thecenter. Nevertheless, many of the same techniques canbe applied to study the physics behind this wide rangeof systems, since the most important ingredients (orbitaland gas dynamics) are common to all of them.Arguably the most convenient experimental test-bedsfor disk dynamics are numerical calculations. The hydro-dynamical equations governing gas dynamics and mag-netohydrodynamics (MHD) have been integrated nu-merically using many approaches. In astrophysics, themost commonly employed techniques are particle-basedmethods like smoothed particle hydrodynamics (SPH,Monaghan & Lattanzio 1985; Springel 2005; Rosswog2009; Price 2012), and grid-based high-resolution shock-capturing techniques (Lax 1957; Godunov 1959; van Leer1977; Woodward & Colella 1984; Brio & Wu 1988;Gardiner & Stone 2005; Mignone et al. 2007), whichuse Godunov-type schemes for hydrodynamic evolutionand often employ adaptive mesh refinement (AMR,Fryxell et al. 2000; Teyssier 2002; O’Shea et al. 2004) toresolve large dynamic ranges. Recently, moving-meshtechniques (Springel 2010; Duffell & MacFadyen 2011;Yalinewich et al. 2015) and several new “meshless” tech-niques (Maron et al. 2012; Hopkins 2015; DeBuhr et al.2015; Hopkins & Raives 2016) have emerged as an at-tempt to merge the accuracy of the AMR approach withthe flow adaptivity of SPH. Moving mesh methods havealready enjoyed remarkable success using an adaptiveVoronoi tessellation for the shape of the mesh zones. On Duffellthe other hand, none of these approaches are specificallytailored to the special challenges inherent in disks.Disks are often highly supersonic (Mach number M & Fig. 1.—
DISCO’s numerical grid is shown in 2D, illustrating therotating annular wedges which make up the cylindrical mesh. Thetest problem being run here is the “Cylindrical Kelvin Helmholtz”test of section 3.1.5. of the grid). RODEO can therefore capture very accu-rate details in the vicinity of a planet, since its orbitalmotion is subtracted and the mesh can also be refinedaround the planet.There is also an efficient code called PEnGUIn(Fung et al. 2014). PEnGUIn is particularly efficient,as it has been optimized for use on graphics processors(GPUs). This is ideal for studies of disks, as many disk-related problems require evolving the system over manythousands of orbits, which is often prohibitively expen-sive. Large parameter surveys are seldom undertakenfor these reasons, but PEnGUIn makes such expensiveproblems much more manageable.While these existing methods are reliable and havedemonstrated themselves useful for evolving disks inmany scenarios, it is worthwhile to consider another dis-tinct numerical approach. In this work, the code DISCOis presented, a moving-mesh technique similar to the nu-merical scheme of the moving-mesh codes AREPO andTESS. In the same spirit of FARGO and RODEO, or-bital motion is subtracted, but unlike RODEO the en-tire shear flow can be subtracted, and unlike the originalFARGO scheme, the topology of the numerical mesh isnot restricted (this last point is particularly important ifone does not wish to excise the inner disk). The com-putational domain is decomposed into zones which arecylindrical wedges (Figure 1). These zones are given anazimuthal velocity and are allowed to shear past one an-other smoothly. This azimuthal velocity can be chosen tohave any value, with the choice of the local fluid velocityresulting in an azimuthally Lagrangian scheme. Becausethe mesh moves with the flow instead of passing the fluidfrom one zone to the next, advection errors resulting fromorbital motion are significantly reduced and otherwisesubtle features can be captured accurately while the floworbits supersonically.DISCO has already been applied to many chal-lenging problems in astrophysics, including gap open-ing and orbital evolution in protoplanetary disks(Duffell & MacFadyen 2012, 2013; Duffell et al. 2014;ISCO Code 3Duffell 2015b; Duffell & Chiang 2015) and the evolutionof circumbinary disks surrounding supermassive blackhole binaries (Farris et al. 2014, 2015b,a; D’Orazio et al.2016). In this work, the numerical technique is describedin detail (Section 2). This includes the recent addition ofa constrained transport method for solving the equationsof magnetohydrodynamics (Section 2.10). In Section 3,a series of numerical code tests is presented, to demon-strate the convergence and practicality of the code. Re-sults are summarized in Section 4. NUMERICAL METHOD
Field Equations
DISCO is capable of evolving arbitrary hyperbolic par-tial differential equations in conservation-law form. Itssimplest mode solves Euler’s equations: ∂ t ( ρ ) + ∇ · ( ρ~v ) = 0 ∂ t ( ρ~v ) + ∇ · ( ρ~v~v + P I ↔ ) = 0 ∂ t ( 12 ρv + ǫ ) + ∇ · (( 12 ρv + ǫ + P ) ~v ) = 0 (1)where ρ is density, ~v is velocity, P is pressure, and ǫ isinternal energy density. An adiabatic equation of stateis typically employed: P = ( γ − ǫ, (2)where γ is the adiabatic index. Additional terms suchas viscosity, gravity, and magnetic fields will be includedin later subsections. For now, the numerical formula-tion will be expressed in terms of these “bare” equations.Other forms, such as the special and general relativisticversions of these equations, will not be discussed here,but will be addressed in a future work.Because angular momentum conservation is so impor-tant to the orbital dynamics, the momentum conserva-tion law is evaluated in terms of the vertical, radial, andangular momentum: ∂ t ( r ρω ) + ∇ · ( r ρω~v + P ˆ φ ) = 0 ∂ t ( ρv r ) + ∇ · ( ρv r ~v + P ˆ r ) = ρω r + P/r∂ t ( ρv z ) + ∇ · ( ρv z ~v + P ˆ z ) = 0 (3)where r is the cylindrical radius, ω is the angular fluidvelocity, and v r and v z are the radial and vertical com-ponents of velocity. The source terms on the right-handside of the equation for radial momentum come from theevaluation of the tensor divergence in cylindrical coordi-nates. This reflects the fact that “radial momentum” isnot a conserved quantity. These terms can be interpretedas a centrifugal force, and an azimuthal pressure-balanceterm that comes out of the non-zero divergence of the ˆ r vector.Finally, before discretizing these equations, the energyequation is modified to improve accuracy. Define thequantity ˜ ω ≡ ω − Ω E ( r ) , (4)where Ω E ( r ) is some differentiable function of radius. Itshould be made clear that Ω E ( r ) has nothing to do withthe mesh motion described in the later sections; it is de-fined in order to subtract off part of the kinetic energyfrom the field equations. For typical applications, Ω E ( r ) Fig. 2.—
Diagram of two adjacent cells in DISCO’s mesh. Theblue shaded region represents a “face” in DISCO; it is defined asthe region of overlap between the surfaces of two neighboring zones. will be the Keplerian orbital velocity, but in principleit can be any analytically known differentiable functionof radius, including Ω E = 0. Nothing in the DISCOalgorithm depends sensitively on the choice of Ω E ( r ),but in practice subtracting off this large kinetic compo-nent of the energy can yield a substantial improvementto the stability and reliability of the numerical scheme,especially for high Mach number flows. Also, note thatΩ E ( r ) will not be subtracted from the angular momen-tum, so that Coriolis terms do not appear in the momen-tum equations.Subtracting Ω E ( r ) from ω in the energy equation yieldsthe following evolution equation (where ˜ v is the velocitywith Ω E subtracted): ∂ t ( 12 ρ ˜ v + ǫ )+ ∇ · (( 12 ρ ˜ v + ǫ + P ) ~v )= rρv r (Ω E ( r ) − r d Ω E dr ˜ ω ) (5)If Ω E ( r ) = 0 is chosen, energy is not explicitly con-served in this formulation, due to the presence of thissource term.To summarize, the field equations can be expressed inthe conservation-law form: ∂ t u + ∇ · ~F = S. (6)The conservation laws can be expressed in terms of thefive primitive variables: W = { ρ, ω, v r , v z , P } . (7)Now, the evolution equations can be compactly summa-rized by writing down the five conserved variables, u = { ρ, r ρω, ρv r , ρv z , ρ ˜ v + ǫ } , (8)the five corresponding fluxes, ~F = { ρ~v , r ρω~v + P ˆ φ , ρv r ~v + P ˆ r ,ρv z ~v + P ˆ z , ( 12 ρ ˜ v + ǫ + P ) ~v } , (9) Duffelland the five source terms, S = { , , ρω r + P/r, ,rρv r (Ω E ( r ) − r Ω ′ E ( r )˜ ω ) } . (10)Now that these have been specified, the evolution ofthe system can be described in terms of the generic ex-pression (6). Mesh Construction Algorithm
Equation (6) will be discretized in the following sub-section (2.3). First, it is necessary to describe DISCO’smesh, and how it is constructed.Similar to the formulation of the TESS code(Duffell & MacFadyen 2011), the numerical scheme canbe completely specified, given the volumes of the zones,logical information specifying which zones are neigh-bors, and the areas of the “faces” connecting neighboringzones. The zones are annular wedges with extents givenin cylindrical coordinates by ∆ r , ∆ φ and ∆ z (Figure2). The faces with ˆ φ normal are the “front” and “back”of these zones (yellow shaded area of Figure 2). Thefaces with ˆ r or ˆ z normal are defined as the overlap of theboundary of two neighboring zones (blue shaded area ofFigure 2). This means that zones can have more thantwo radial or vertical faces (on average, zones typicallyhave four of each).Given neighboring annuli at radii r j and r j +1 , firstthe zone at each of these radii intersecting the radial ray φ = 0 are found. These two zones are guaranteed to sharea face. The geometry of the shared face is identified, andthen the next face is found by advancing whichever of thetwo zones has a smaller φ i +1 / associated with their frontface. Again, this new pair of zones must share a face.This procedure is repeated for N φ ( j ) + N φ ( j + 1) faces,where N φ ( j ) is the number of zones in a given annuluslabeled by the index j . Each step, a face is identified asthe intersection of the boundary of the pair of zones.Faces normal to ˆ z are constructed in an analogous way. Integral Form / Mesh Motion
DISCO is a finite-volume method. To discretize thesystem, (6) is integrated over the volume of a computa-tional zone, using Gauss’ law on the flux term: Z ∂ t udV + I ~dA · ~F = Z SdV. (11)Now, at first, consider the case of no mesh motion, forwhich Z ( ∂ t u ) dV = ∂ t Z udV, (12)and define M ni to be the amount of a given conservedquantity in zone i at timestep n: M ni ≡ Z udV (13)After performing an integral in time, it is then straight-forward to write an evolution equation for the { M ni } : M n +1 i = M ni − ∆ t X face f ~dA f · ~F f − S i ∆ t ∆ V. (14) where the sum is over faces bounding zone i . So far,no approximations have been made, so long as ~F f is in-terpreted as the time-averaged and face-averaged flux,and S i is interpreted as the volume-averaged and time-averaged source term. Exact geometry is employed, sothat for example ∆ V is the exact volume of a cylindricalwedge: ∆ V = ∆ φ ∆ z ( 12 r − r − ) , (15)where ∆ φ and ∆ z are the azimuthal and vertical extentof the zone, and r + and r − are the outer and inner radii.Some of the conserved quantities and source terms re-quire the coordinate r. In this case, the radius is chosento be the moment arm of the zone: r moment = r
12 ( r + r − ) (16)If the orbital motion of the grid is switched on, the inte-gral form of these equations departs from (14), since thecontrol volumes and their associated faces move throughspace. In this case, the generalization of Equation (12)is given by the Reynolds transport theorem as Z ( ∂ t u ) dV = ∂ t (cid:18)Z udV (cid:19) − I u ~w · ~dA, (17)where ~w is the velocity of the boundary of the zone. Thisresults in the following modification to Equation (14): M n +1 i = M ni − ∆ t X face f ~dA f · ( ~F f − ~w f u f ) − S i ∆ t ∆ V. (18)In other words, the time evolution differs from a stan-dard Godunov-type method by making the substitution ~F → ~F − ~wu , where ~w is the velocity of the face uponwhich the flux is evaluated, and now u on the face needsto be determined in addition to F . For DISCO’s cylin-drical mesh, the velocity ~w is always normal to the face.The face velocities are zero for all radially and verticallyoriented faces; all mesh motion is azimuthal. The az-imuthal velocity can be chosen in many ways. If desired,each azimuthal face can be moved independently, accord-ing to the average velocity of its adjacent zones. Alter-natively, some global analytical formula for ~w ( r ) can beprescribed, if one does not wish for the zones to have toomuch independence. The mesh velocity could also be setto w = r Ω E ( r ), where Ω E ( r ) is the analytical functiondescribed in the previous subsection. Again, it is notnecessary to force these two functions to be equal.So far, no numerical approximations have been made;equation (18) is merely an integral form of (6). There-fore, the numerical approximations are housed in the es-timation of ~F f , the time-averaged and area-averaged fluxthrough the face (as well as u f and S i ). These numericalapproximations are detailed in the following two subsec-tions. Riemann Solver
Equation (18) requires a numerical estimate for F f ,the time-averaged and area-averaged flux through faceISCO Code 5 (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:2) (cid:1) (cid:1) (cid:3) (cid:1) (cid:4) (cid:2) (cid:5) (cid:6) Fig. 3.—
Schematic diagram of the Riemann problem on a mov-ing mesh. The face velocity is traced out by the red dashed curve x = wt ; this is the characteristic on which the solution is evaluated. f . The standard method for calculating such a flux inthe Godunov method is to use a Riemann solver.A Riemann solver takes as input a left and right state { u L } , { u R } and returns as output some estimate of thesolution to the shock-tube problem given by piecewise-constant initial conditions: u ( x, t = 0) = (cid:26) u L x < u R x > x = 0: F ∗ = F ( x = 0 , t ) (20)In the moving-mesh case, a different output is desired.If the interface moves with velocity w , then the fluxshould be evaluated along the characteristic x = wt : F ∗ = F ( x = wt, t ) (21) u ∗ = u ( x = wt, t ) (22)(see Figure 3). Evaluating the Riemann solution along agiven characteristic ( x = wt instead of x = 0) is straight-forward.DISCO employs several different approximate Rie-mann solvers. HLLE and HLLC (Toro 2013) areavailable for all flows, though HLLC is necessaryto preserve contact discontinuities to high precision(Duffell & MacFadyen 2011). For MHD flows, an HLLDsolver is implemented in DISCO, based on the solverof Miyoshi & Kusano (2005). Analogous to HLLC, theHLLD solver is necessary for preserving magnetic discon-tinuities to high precision (advection of a field loop canbe solved to machine precision, but only if the HLLDsolver is employed). Piecewise Linear Reconstruction
In order to achieve second-order accuracy in space,primitive variables must be extrapolated from zone cen-ters to faces to produce input to the Riemann solver (Fig-ure 4):
Fig. 4.—
Extrapolation to the face adjoining two cells. Themisalignment of the cells means that the gradients in the azimuthaldimension are needed to extrapolate primitive variables to radiallyand vertically oriented faces. W L = W i + ( ~x f − ~x i ) · ( ~ ∇ W ) i (23) W R = W j + ( ~x f − ~x j ) · ( ~ ∇ W ) j (24)where the ~ ∇ W are slopes which are estimated at the zonecenters. In calculating these slopes, care must be takenas the primitive variables cannot be assumed to representdifferentiable functions. Therefore, after estimating gra-dients of the primitive variables, a slope-limiter is appliedto ensure stability in non-smooth regions of the flow.First, the azimuthal gradients are calculated. This isdone by calculating left, right, and centered gradients inthe zone: S L = ( W i − W i − ) / ( r ∆ φ L ) (25) S R = ( W i +1 − W i ) / ( r ∆ φ R ) (26) S C = ( W i +1 − W i − ) / ( r (∆ φ L + ∆ φ R )) (27)where ∆ φ L = (∆ φ i + ∆ φ i − ) /
2, and ∆ φ R = (∆ φ i +∆ φ i +1 ) /
2. The slope-limited azimuthal gradient is thengiven by ∇ φ W = minmod( θ plm S L , θ plm S R , S C ) , (28)where θ plm is a slope-limiting parameter 1 < θ plm < x, y, z ) = ( min( x, y, z ) x, y, z > x, y, z ) x, y, z <
00 otherwise (29)Next, the radial and vertical gradients are calculated.For brevity, the formulas for the vertical gradients areomitted, as they are similar to the formulas for the radialgradients.First, the radial gradient at each radially-oriented faceis estimated, using the extrapolated values from the az-imuthal gradient: Duffell W if ≡ W i + r i ∆ φ if ( ∇ φ W ) i (30) h∇ r W i face ij = W if − W jf r i − r j , (31)where ∆ φ if is the angular separation between the centerof zone i and the center of face f . The zone-centeredgradient is then estimated by performing an average overfaces, weighted by face area: h∇ r W i zone i = P j dA j h∇ r W i face ij P j dA j (32)This provides a radial gradient in each zone, whichlike the azimuthal gradient must be processed through aslope limiter for stability. This slope limiter is essentiallythe same as the one used in the azimuthal direction, butincluding all neighbors. The “centered” slope has alreadybeen calculated above; the final slope used is given by ∇ r W = minmod( h∇ r W i zone i , n θ plm h∇ r W i face ij o ) . (33)Formulas (30)-(33) are then repeated for vertical gra-dients so that the gradient ~ ∇ W is fully determined in thezone. For most problems, the slope-limiting parameteris chosen to be θ plm = 1 . Time Evolution
Equation (18) specifies how to advance from timestep n to timestep n + 1, given the time-averaged values of F and u on the face. This is given by the Riemannsolver, which takes as input a left and right state, { W L } and { W R } . This left and right state are found by ex-trapolating from zone centers to face centers, using theslope-limited gradients given in section 2.5.At this point, the system is completely specified, butthe time-evolution operator, as expressed in equation(18), is only first-order in time. It has the following form: M n +1 i = M ni + ∆ tL i ( { state n } ) , (34)where L is a time-evolution operator depending on thestate of the system at timestep n . To increase the orderof accuracy of the code, a method-of-lines technique isemployed, introducing the intermediate state M (1) : M (1) i = M ni + ∆ tL i ( { state n } ) , (35) M n +1 i = 12 ( M (1) i + M ni ) + 12 ∆ tL i ( { state (1) } ) . (36)This constitutes a second-order timestep which is con-sistent with a total variation diminishing scheme. Thetimestep ∆ t is Courant-limited; that is,∆ t < min(∆ t cross i ) , (37)where ∆ t cross is the shortest signal-crossing time of azone: Fig. 5.—
Parallelization is accomplished by domain decompo-sition in the radial and vertical dimensions. Two cut planes aredisplayed in 3D showing which processor each zone belongs to.Each color represents a single processor (in this example, the workis divided among 25 processors). ∆ t cross = min (cid:18) ∆ rc s + | v r | , ∆ zc s + | v z | , r ∆ φc s + | v φ − w | (cid:19) . (38)Note that moving the zones to cancel a supersonic or-bital velocity ( w ∼ v φ ≫ c s ) can increase the allowedtimestep by orders of magnitude. In the MHD case thesound speed in the above formula is replaced by the speedof fast magnetosonic waves.The courant condition (37) is satisfied at each time-step by setting ∆ t = C CFL ∆ t cross , where C CFL < . . Parallelization
DISCO achieves efficient parallelization by subdividingthe computational domain into annuli. Define N Gr and N Gz to be the global radial and vertical dimensions of thecomputational grid. Define the indices n r and n z whichlabel the radial and vertical grid: 0 < n r < N Gr and0 < n z < N Gz . The number of zones in the azimuthaldimension can vary with r and z : N φ = N φ ( n r , n z ).Typically N φ is chosen at each radius so that the zoneshave a nearly 1 : 1 aspect ratio. The domain is thensubdivided in r and z so that each processor has a localresolution N Lr , N Lz (Figure 5). Boundary data is shippedin the vertical and radial direction every timestep.It is also possible to subdivide the domain in azimuth,but such a subdivision adds significant complexity to themethod. The official version of DISCO therefore onlyperforms parallel subdivision of the domain radially andvertically (in other words, N Gφ = N Lφ = N φ ). Perfor-mance and scaling of DISCO on thousands of CPUs istested in section 3.5. Disk-Satellite Interactions
Disks usually orbit around a central point mass, mean-ing gravitational source terms must enter the evolutionequations. Additionally, there may be orbiting satellitesexerting a gravitational influence on the disk. These in-fluences are accounted for by adding the following sourceterms to the energy and momentum equations: S Momentumgrav = ρ~g, (39)ISCO Code 7 S Energygrav = ρ ( ~g · ~v ) (40)where ~g is the gravitational acceleration due to all pointmasses, labeled by p : ~g = X p −∇ Φ p (41)Φ p = − Gm p ( | ~x − ~x p | n + ǫ np ) /n (42)Here, n and ǫ p are optional smoothing parameters.Typically in 2D calculations, n = 2 and ǫ p = 0 . h , where h is a scale height. This is to mimic the averaging of thegravitational force over the disk scale height (if there isa central body, typically ǫ p = 0 for that body). In thecode, all point masses (e.g. stars, planets, black holes)are simply called “planets”, and treated identically in analgorithmic sense. The point masses can be given anyprescribed motion ~x p ( t ), or be moved due to the gravi-tational influence of the gas.Accretion onto these bodies is also possible using anadditional source term in the continuity equation. Thishas been employed in studies of binary-disk interactions(Farris et al. 2014, 2015b,a; D’Orazio et al. 2015), but isnot in DISCO’s public version as there are many possiblechoices for such a term. Viscosity
A Navier-Stokes viscosity can be represented as asource term in the momentum equation, ~S = ∇ · σ ↔ , (43)and a source term in the energy equation, S = ( ∇ · σ ↔ ) · ˜ v + σ ij ∇ i v j , (44)where the viscous stress tensor σ ↔ will be defined below.The first term in the energy equation can be interpretedas work done by viscous forces (inner product of forcewith velocity, F · v ) and the second term expresses viscousheating. Both of these source terms can be re-expressedas a viscous flux: F Momentumvisc = − σ ↔ , (45) F Energyvisc = − σ ↔ · ˜ v, (46)This is possible because viscosity is an internal body-force in the gas, and therefore conserves total momentumand energy. In the case that Ω ′ E ( r ) = 0, the energyequation has a source term: S Energy visc = σ rφ r Ω ′ E ( r ) . (47)In cylindrical coordinates, the tensor divergence gen-erates a more complicated expression for the fluxes, in-cluding a source term for radial momentum. The viscousstress tensor σ ↔ is proportional to the velocity gradients: σ ij = νρ (( ∇ i v j + ∇ j v i ) + ηδ ij ∇ · v ) (48) η is a dimensionless order-unity constant which sum-marizes the relationship between bulk and shear viscos-ity. For most orbital flows, the choice of η is unimportant (cid:1) (cid:1) (cid:2) (cid:3) (cid:4) (cid:5) (cid:2) (cid:6) (cid:7) (cid:3) (cid:1) (cid:8) Fig. 6.—
Schematic diagram of Faraday’s law in integral form;the flow of magnetic field lines across an edge is equal to the lineintegrated electric field along the edge. (in fact, for most problems it will only be the value of σ rφ that matters). The full derivation appears in theappendix; the final values of the viscous flux and sourceterms are ~F visc = − νρ { , r ~ ∇ ω + 2 v r ˆ φ, ~ ∇ v r − ω ˆ φ, ~ ∇ v z ,v r ~ ∇ v r + r ˜ ω ~ ∇ ω + v z ~ ∇ v z − rω ˜ ω } , (49) S visc = { , , − νρv r /r , , − ρν ( ∇ φ v r + r ∇ r ω ) r Ω ′ E ( r ) } . (50)where ˜ ω = ω − Ω E ( r ). In section 3.2, several test prob-lems will be presented to empirically check that all ofthese terms are correct.Note that the formula (49) is expressed in terms ofgradients of the primitive variables. The viscous fluxis evaluated separately from the Riemann solver. Theprimitive variables W are extrapolated to the face to at-tain the quantities W L and W R , and these quantities andtheir slope-limited gradients ~ ∇ W are averaged betweenthe left and right state, before using them to evaluatethe viscous fluxes (49). Magnetohydrodynamics
Similar to Euler’s equations, the field equations ofMHD can also be expressed in conservation-law form: ∂ t ( ρ ) + ∇ · ( ρ~v ) = 0 ∂ t ( ρ~v ) + ∇ · ( ρ~v~v + ( P + 12 B ) I ↔ − ~B ~B ) = 0 ∂ t ( 12 ρv + ǫ + 12 B )+ ∇ · (( 12 ρv + ǫ + P + B ) ~v − ( v · B ) ~B ) = 0 ∂ t ( ~B ) + ∇ · ( ~v ~B − ~B~v ) = 0 (51)Rewriting the equations in a cylindrical basis, it is pos-sible to express these as additions to the standard hydrovariables: u = u hydro + u mhd , (52)where the u hydro are given by (8), and u mhd = { , , , , B , B r , B φ /r, B z } . (53) DuffellThe MHD fluxes are similarly summarized as ~F mhd = { , r ( 12 B ˆ φ − B φ ~B ) , B ˆ r − B r ~B, B ˆ z − B z ~B, B ~v − v · B ~B,B r ~v − v r ~B, ( B φ ~v − v φ ~B ) /r, B z ~v − v z ~B } . (54)and there is a new source term for radial momentum,given by magnetic tension, or “hoop stress”: S mhd = { , , B φ /r, , , , , } . (55)MHD is a challenging set of field equations to inte-grate, because of subtle behaviors relating to the diver-gence constraint, ∇ · B = 0 (e.g. Brackbill & Barnes1980; Evans & Hawley 1988; Dai & Woodward 1998;Ryu et al. 1998; Balsara & Spicer 1999; T´oth 2000).One way of framing the ∇ · B problem is in termsof the following thought experiment: consider a two-dimensional cartesian domain with a field loop advectingwith a uniform velocity in the x direction. While B y isadvected, there is no x-directed flux for the parallel com-ponent B x . Instead, the advection of B x is accounted forby the rotational flux in the y direction. ∂ t B x − ∂ y ( v x B y ) = 0 (56)but the divergence constraint implies ∂ y B y = − ∂ x B x (57)so if ∇ · B is guaranteed to be zero by the numerical sten-cil which updates the magnetic field, then the rotationalflux in y is identical to an advective flux in the x direc-tion, so it makes no difference that B x is not explicitlyadvected in this sense.However, if one cannot guarantee that ∇ · B = 0 thensomething as simple as advecting a field loop can gowrong. B x and B y are numerically updated in a funda-mentally different way, and this causes the loop to even-tually destabilize. This becomes more of a problem thelarger ∇ · B is allowed to grow.Many techniques have been developed to stably evolvethis system, the most successful of which is constrainedtransport (CT, Evans & Hawley 1988). However, CT isoften difficult to implement on complicated grids, in par-ticular moving meshes such as DISCO’s (though notably,Mocz et al. (2014) employed CT on a 2D Voronoi mesh,demonstrating that CT is possible even with complex ge-ometries). As a result, several methods have been devel-oped which have less dependence on the mesh employed.One popular such technique is to modify the evolutionequations so that the divergence constraint propagatesand diffuses (e.g. Dedner et al. 2002).Such divergence-cleaning techniques are easy to imple-ment, but it is difficult to test their effectiveness, becausethey do not guarantee machine-precision zero divergencefor any stencil. Even a small divergence error can causeinaccurate physics on long timescales.Another issue with the formulation of Dedner et al.(2002) is that it introduces an additional wavespeed intothe system. In the standard formulation, this wavespeedmust be the fastest velocity in the system, in order thatthese waves can keep up with divergence errors quickly enough to correct them. Unfortunately for the movingmesh technique, this eliminates the time-step advantage,because this wave moves quickly with respect to the grid.Alternate Galilean-invariant formulations of theseequations are possible (e.g. Powell et al. 1999), but thesealways require source terms which have derivatives,which undermines many of the advantages of the finite-volume formulation (one principal advantage of finite vol-ume methods is that they evolve the integral form of theequations, and as a result do not necessitate smooth solu-tions). In short, it may be more advantageous to preventdivergence errors from appearing in the first place.Another method increasingly employed is the vec-tor potential formulation (Del Zanna et al. 2003;Etienne et al. 2010). This, of course, has the advantagethat divergence errors are never introduced, as B isdefined as a curl. It is also a preferred method fornontrivial meshes, which can complicate implementa-tions of CT. Some vector potential formulations havebeen shown to be functionally equivalent to CT onuniform meshes (Etienne et al. 2010; Helzel et al. 2011;Helzel et al. 2013).Unfortunately, the vector potential formulation mightalso require a numerical derivative, in the operation B = ∇ × A , and therefore some formulations might as-sume that the vector potential is differentiable (for ex-ample, this is generally the case if the vector potential iscell-centered). Another disadvantage is that introducinga vector potential also introduces gauge modes. Eitherthese gauge modes are static, and accumulate on the grid,or they propagate at some velocity which must be intro-duced, and they can also limit the code’s time-step, sim-ilar to Dedner’s method. The time-step advantage is bigenough for the moving mesh that it is worth maintaining,if possible. Of course, none of these disadvantages canbe truly devastating for vector potential formulations, asany constrained transport scheme could be re-expressedas a vector potential scheme, by evolving ~A on each edgeof the mesh using the same electric fields, and calculat-ing the magnetic flux through each face by integrating A around a closed loop. In this case, the truncation errorcan be identical to a CT scheme.Constrained transport techniques find ways of process-ing the MHD fluxes so that they do not introduce anydivergence errors. In a sense, the idea is somewhat anal-ogous to conservative formulations, which evolve the sys-tem in such a way as to avoid conservation-law violations.The CT formulation of Evans & Hawley (1988), themost commonly employed CT scheme, makes use of thenatural topology of the MHD equations. Faraday’s lawcan be expressed as a conservation law, and it is straight-forward to define conserved quantities in volumes, andfluxes through faces. On the other hand, ~B is a conservedflux . In other words, Faraday’s law is more naturally ex-pressed in a lower-dimensional form, by integrating itover the area of a face: ∂ t Φ = − I ~E · ~dl. (58)This is a lower-dimensional analog to the finite-volumeconservation law (equation 11). The magnetic flux Φ isthe analog of mass, and the electric field E is the analogof mass flux (the direction of the flow is perpendicular toISCO Code 9 (cid:1) (cid:1) (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:4) (cid:1) (cid:5) (cid:1) (cid:2) (cid:1) (cid:1) (cid:2) (cid:2) (cid:1) (cid:2) (cid:3) (cid:1) (cid:2) (cid:4) (cid:1) (cid:3) (cid:5) (cid:1) (cid:3) (cid:6) (cid:1) (cid:3) (cid:7) (cid:1) (cid:3) (cid:8) Fig. 7.—
Upper Panel: Zone-specific faces are identified. Eachzone can be uniquely associated with five faces in three dimensions(or three faces in two dimensions). Lower Panel: Zone-specificedges are identified. Radial and vertical edges are associated withthe interior a given zone. Some of these edges are redundant, asthey are duplicated in other zones. This is accounted for by av-eraging the electric field between associated edges on neighboringzones. the electric field, but the component of E along the edgeprovides the magnitude of the flux).Given a computed value of ~E = − ~v × ~B , one can gen-eralize the integral equations derived in section 2.3, toarrive at the “finite-area” form of (58):Φ n +1 i = Φ ni − ∆ t X edge e ~dl e · ~E e . (59)The numerical scheme can be interpreted as a lower-dimensional analog of standard finite-volume methods,with Φ playing the role of mass, and E being the ana-log of mass flux (see Figure 6). Because the fluxes Φ arestored on faces, this requires some memory of the numer-ical mesh to persist from one time-step to the next. Forconvenience, each flux Φ is stored on one of the zoneswhich houses the face associated with this flux (see Fig-ure 7). The 3D CT scheme used in DISCO is based onthe method of Evans & Hawley (1988), though the elec-tric fields are not computed in exactly the same way. Face-Centered Fluxes to Cell-Centered Fields
The complete description of the time-evolution for themagnetic fluxes will appear in section 2.10.2. First, it isimportant to recognize that it will not be possible to cal-culate the magnetic field directly from the magnetic flux.This would require the operation h B i = Φ /dA , where dA is the area of the face. Unfortunately, this area can bearbitrarily small, as faces can disappear or reappear dur-ing changes in the topology of the mesh. Calculating themagnetic field by dividing a small flux by a small areacan result in arbitrarily large errors; small errors in theflux would lead to arbitrarily large errors in the field.Therefore, in DISCO’s numerical scheme the operation h B i = Φ /dA is never explicitly performed. Instead, thefluxes on each face are used to determine a cell-centeredaverage magnetic field, and this magnetic field is usedto update the system. In fact, most of the timestep pro-ceeds as a normal finite-volume scheme, treating the cell-centered magnetic fields as primitive variables; they areinterpolated to the faces along with the other variables,and used as input to Riemann solvers. The only time theface-centered fluxes are needed is in the step just beforeconverting from conserved to primitive variables, wherethe cell-centered magnetic fields are determined based onthe face-centered fluxes. In this sense, one could interpretthe update of the magnetic fields as a predictor-correctorscheme (Helzel et al. 2011).This operation (from Φ face to B zone ) is performed bydetermining the total flux through each surface of thezone; the flux is then assumed to vary linearly throughthe zone, so that h ~B · ˆ n i cell dA c = 12 ( X side 1 Φ + X side 2 Φ) . (60)In other words, the cell-centered components of themagnetic field are h B r i cell = ( P Φ in + P Φ out )( r + + r − )∆ φ ∆ z (61) h B z i cell = ( P Φ top + P Φ bottom )( r + + r − )∆ φ ∆ r (62) h B φ i cell = (Φ front + Φ back )2∆ r ∆ z (63)Once these zone-centered fields are calculated, the con-served variables in the zone are converted to primitivevariables, using the zone-centered magnetic fields justcomputed. The rest of the time-update proceeds usingthese fields as cell-centered primitive variables. Time-Update of the Magnetic Fluxes
In section 2.10.3, it will be explained how the electricfield is calculated on each edge. For now, it is assumedthat this E field is known, and the goal is to derive a ruleto update the magnetic fluxes during a timestep. Thisstarts with Faraday’s law: ∂ t ~B + ∇ × ~E = 0 (64)Following the same process which led to (18), Faraday’slaw is integrated over the surface of a face: ∂ t Φ + I ~E · dl = 0 (65)Now, interpreting ~E e as the time-averaged electric fieldon edge e , the time-update step can be expressed as:0 DuffellΦ n +1 f = Φ nf − ∆ t X edge e ~E e · dl e . (66)This method of evolving face-centered fluxes doesnot conserve volume-integrated flux density, as someuniform-mesh CT techniques do (e.g. T´oth 2000). If theedges move with velocity ~w , the generalization is:Φ n +1 f = Φ nf − ∆ t X edge e ( ~E e + ~w e × ~B e ) · dl e . (67)This is also how mesh motion is accounted for byMocz et al. (2014). Just like equation (18), this is an ex-act expression, suitably interpreted. All of the numericalapproximations will be housed in the calculation of E e and B e , the time-averaged electric and magnetic fields.These will be determined in the following section. Calculation of the Electric Fields
Any choice for the electric fields will preserve the di-vergence constraint, but it is important to make a well-motivated choice in order to provide an accurate and sta-ble approximation to the underlying field equations. Forexample, Balsara & Spicer (1999) showed that simplyaveraging the electric fields on four adjacent faces failsto produce an upwind scheme, and as a result this cangive innacurate solutions. Various remedies are proposedfor this (e.g. Balsara & Spicer 1999; Gardiner & Stone2005). Here, the numerical scheme is not as suscepti-ble to this upwinding issue, because the mesh motiontypically keeps azimuthal faces within the Riemann fan.However, it is still nontrivial to develop a stable schemebecause of the complicated mesh structure. In practice,a stable scheme was found by experimenting with differ-ent averaging procedures. The calculation of the electricfields will involve three steps: first, a definition of zone-specific edge-centered fields, secondly an identification ofthese electric fields with fluxes found in the Riemannsolver step, and finally an averaging process over severalneighboring zones.Indexing of the faces and edges in each zone is shownin Figure 7. To calculate the electric field on these edges,note that this electric field gives the number of magneticfield lines being dragged across each edge per unit time.If one knows, for example, the flow of B φ through a facewith r-directed normal, then one implicitly knows theline integral of E z along a z-directed edge parallel to thatface.These statements can be made more mathematicallyconcrete. Using the geometrical fact that the area ele-ment is given by the wedge product of two line elements dA i = ǫ ijk dx j ∧ dx k , the integrated flux of a magneticfield component B i through a face is given by F ( B i ) · dA = ( v j B i − v i B j ) 12 ǫ jkl dx k ∧ dx l (68)The combination ( v j B i − v i B j ) can be related to thecross product ~v × ~B : F ( B i ) · dA = ǫ jim ( v × B ) m ǫ jkl dx k ∧ dx l (69) Given ~E = − ~v × ~B and the identities of the ǫ pseu-dotensor, one can arrive at the relation: F ( B i ) · dA = E j dx j ∧ dx i . (70)If one integrates this formula over a face, then dividesby area, one attains h F ( B i ) i = R E · dl ∧ dx i dA , (71)where brackets denote area-averaged fluxes.Therefore, the line-integrated electric field can be re-lated to the area-integrated flux of magnetic field, whichhas already been calculated from the Riemann solver. Inother words, one can associate a line-averaged E withan area-averaged flux. For each electric field component,there are two areas to average over, corresponding to thetwo coordinate planes parallel to this edge. The electricfield is set to the mean of these two area-averaged fluxes.These fields are not averaged with neighboring zones yet,as this will be done in the next step. Specifically, for thevertical fields, E z = 12 D ~F B r · ˆ φ E front face − D ~F B φ · ˆ r E inner face (72) E z = 12 D ~F B r · ˆ φ E front face − D ~F B φ · ˆ r E outer face (73) E z = 12 D ~F B r · ˆ φ E back face − D ~F B φ · ˆ r E inner face (74) E z = 12 D ~F B r · ˆ φ E back face − D ~F B φ · ˆ r E outer face (75)where “front” implies the face with greater value of φ .For the radial fields, E r = 12 D ~F B φ · ˆ z E bottom face − D ~F B z · ˆ φ E front face (76) E r = 12 D ~F B φ · ˆ z E top face − D ~F B z · ˆ φ E front face (77) E r = 12 D ~F B φ · ˆ z E bottom face − D ~F B z · ˆ φ E back face (78) E r = 12 D ~F B φ · ˆ z E top face − D ~F B z · ˆ φ E back face . (79)Azimuthal electric fields will be treated separately, butwill have the same essential forms: E φ = 12 D ~F B z · ˆ r E − D ~F B r · ˆ z E (80)This choice of electric field is designed for consistencyof the numerical scheme, since the flux of field lines isconsistent with the Godunov fluxes calculated from theRiemann solver.After calculating these zone-specific edge-centeredelectric fields, an averaging process is performed betweenzones in order to acquire a well-posed electric field oneach individual edge in the mesh. Self-Consistent Averaged Electric Fields
Once zone-specific electric fields have been specified viaequations (73) - (79), an averaging process is performed,ISCO Code 11
Zone i+1Zone i (cid:1) (cid:2) (cid:1) (cid:1) (cid:2) (cid:2) (cid:1) (cid:2)(cid:3)(cid:4)(cid:5)(cid:6)(cid:7)(cid:8)
Zone j+1Zone j (cid:1) (cid:2) (cid:3)
Fig. 8.—
Neighboring zones have inconsistent electric fields.These electric fields are made consistent by replacing each withthe average of the two. In the upper panel, edges are identifiedbetween zones which share a φ -normal face (equations 81-84 and89-92). In the lower panel, the interpolation procedure is shown(corresponding to equation 86), to identify edges in zones on neigh-boring annuli. in order to ensure that electric fields on adjacent facesare consistent with one another. For example, as Figure8 shows, E z on zone i should be identical to E z on zone i + 1. Before the averaging process, these fields are notgenerally consistent with one another. After averaging,the electric field is self-consistent and the MHD updatestencil effectively spans a larger number of zones.As Figure 8 demonstrates, E z and E z of zone i mustbe compatible with E z and E z of zone i + 1. This isaccomplished via the substitution: E z i →
12 ( E z i + E z i +1 ) (81) E z i →
12 ( E z i + E z i +1 ) (82) E z i →
12 ( E z i + E z i − ) (83) E z i →
12 ( E z i + E z i − ) (84)Similarly, E z of zone j must be compatible with theinterpolated electric field of zone j + 1: (cid:1) (cid:2)(cid:3) (cid:4) (cid:5) (cid:1) Fig. 9.—
Azimuthal electric fields are found by averaging overfour adjacent faces. The resultant electric field is calculated inequation (98). E z j →
12 ( E z j + E interp z j − ) (85) E z j →
12 ( E z j + E interp z j +1 ) (86) E z j →
12 ( E z j + E interp z j − ) (87) E z j →
12 ( E z j + E interp z j +1 ) (88)where E interp z j +1 for example is a weighted average of E z and E z interpolated to the position of E z j or E z j (Figure 8).Similarly, for the radial electric fields: E r i →
12 ( E r i + E r i +1 ) (89) E r i →
12 ( E r i + E r i +1 ) (90) E r i →
12 ( E r i + E r i − ) (91) E r i →
12 ( E r i + E r i − ) (92) E r k →
12 ( E r k + E interp r k − ) (93) E r k →
12 ( E r k + E interp r k +1 ) (94) E r k →
12 ( E r k + E interp r k − ) (95) E r k →
12 ( E r k + E interp r k +1 ) (96) E φ is not averaged with neighbors in this way. Rather,the φ component of the electric field is not defined zone-wise, but on the intersection between four faces (Figure9). In this case, E φ = 14 ( h F ( B z ) · ˆ r i U + h F ( B z ) · ˆ r i D (97) − h F ( B r ) · ˆ z i L − h F ( B r ) · ˆ z i R )Note that on average each zone is in contact with 16azimuthal edge segments, and therefore E φ is effectivelyaveraged over more than four faces.2 Duffell (cid:1) (cid:1) (cid:1) (cid:2) (cid:1) (cid:3) (cid:1) (cid:1) (cid:1) (cid:2) (cid:1) (cid:3) Fig. 10.—
Schematic diagram of a topology flip in the mesh.After a topology flip, magnetic fluxes Φ , Φ and Φ are adjustedfor the new topology using equations (99 - 101) Topology Changes
Because the zones shear past one another, faces candisappear or emerge during the course of a timestep (Fig-ure 10). Ordinarily, this might be a problem for the nu-merical method; when a face disappears, where does itsflux go? When a face emerges, what sets its flux?The way this scheme is designed, the numerical solu-tion is not very sensitive to the answers to these ques-tions, since the area of the face in question is very small inthese circumstances, and the cell-centered magnetic fieldis an average over faces, weighted by face area. Nonethe-less, some flux must be specified. Since the face is chang-ing topology, its bounding edges must be very close to-gether and therefore the electric and magnetic fields mustbe very similar (this is guaranteed by the interpolationand averaging process). Therefore, most of the time-dependence of the flux must come from the difference inthe velocity ~w at each edge:∆Φ ≈ ( w L − w R ) Bdldt. (98)In other words, during a topology flip the change in fluxin the cell is assumed to be entirely due to the amountof flux overtaken by the edge. In this case, during thetime update, there needs to be an adjustment in order toaccount for the fact that in the first part of the timestep,the face was giving up its flux to one zone, and in thelast part of the timestep, it was taking flux from anotherzone.Note that because | w L − w R | dldt > dA during a time-step where the topology flips, Equation (98) guaranteesa sign change in Φ during a topology change; the point atwhich Φ goes to zero signifies the topology flip. There-fore, in Figure 10 all flux lost before the sign changeshould be given to face 1 and all flux accumulated afterthis point should be taken from face 2.Algorithmically, all of this is simple to account for.First, a normal timestep is taken. Then, after thetimestep, if a face “flips”, i.e. a face disappears from onezone and moves onto another zone, then the followingcorrections must be made to this face and the adjacentfaces (see Figure 10):Φ new0 = − Φ (99)Φ new1 = Φ + Φ (100)Φ new2 = Φ + Φ (101) This adjustment also guarantees that the magnetic di-vergence remains zero for the appropriate stencil. Thismethod for tracking topology changes is different fromMocz et al. (2014), where flux is redistributed equallybetween neighboring faces. TEST PROBLEMS
Given the uniqueness of this numerical scheme, a widerange of tests is performed. Many of these tests are“sanity checks”, ensuring that all terms are correctly ac-counted for in the code. Accuracy and convergence arealso very important, and therefore several convergencetests are performed. Additionally, several challengingtests relevant to astrophysics are performed, in orderto demonstrate the code’s robustness and usefulness forstudying nontrivial astrophysical flows, including tests ofdisk-planet interactions, and MHD turbulence driven bythe magnetorotational instability in an accretion flow.Nearly all tests use an adiabatic index of γ = 5 / γ = 1 .
4, and the supersonic Keplerian spreading test,which uses a nearly isothermal adiabatic index γ = 1 . θ plm = 1 . . . Hydrodynamics
Cylindrical Shock Tube
This test demonstrates DISCO’s ability to capture ra-dially propagating shocks. For characteristics movingradially, DISCO performs as a robust high-resolutionshock-capturing code, even when the mesh motion is notutilized.The domain extends from 0 < r <
1, with uniformradial resolution. An adiabatic index of γ = 5 / r < . ρ = 1 . P = 1 .
0. Otherwise, ρ = 0 .
125 and P = 0 .
1. Initially the fluid is not moving: v r = ω = 0. This calculation is run until a time t = 0 . Cylindrical Isentropic Pulse
The isentropic pulse is a simple nonlinear test for con-vergence of any code. Convergence is checked by measur-ing entropy conservation. The initial setup is as follows:The domain extends from 0 < r < .
5, with uniformradial resolution. Density is given by ρ = 1 + 3 e − (80 r ) . (102)Pressure is chosen to be isentropic: P = ρ γ (103)and the fluid is initially stationary, v r = ω = 0. Thepulse explodes outward (Figure 12), and eventually formsa shock, but before the shock forms the equation P = Kρ γ continues to hold due to entropy conservation ( K =ISCO Code 13 Fig. 11.—
Cylindrical shock tube test at t = 0 .
25. DISCO’s gridwith 128 radial zones is compared with an identical high-resolutioncalculation using a 1D code in cylindrical coordinates. Qualitativeagreement is found; DISCO accurately captures the shocks prop-agating radially, even though mesh motion is not utilized in thistest.
Fig. 12.—
Cylindrical isentropic pulse test at t = 0 . N r = 64 radial zones. The pulse is offset with respect to theorigin to test the convergence of a problem which includes all hydrofluxes. In this test, each zone is moved individually, so that zonesare compressed or expanded with the flow. P/ρ γ evolves as a conserved scalar, as long as the solutionremains smooth).Error is calculated by verifying entropy conservationat a time t = 0 .
1, before the shock has formed: L = R | P/ρ / − . | dV R dV . (104)Fast convergence is found for this problem (Figure 13); D en s i t y -7 -6 -5 -4 -3
100 1000N -2 L E rr o r ResolutionCentered on originOffset, with Mesh Motion
Fig. 13.—
Convergence of the cylindrical isentropic pulse test.The upper panel compares the run with N r = 64 with a high-resolution calculation with a 1D code in cylindrical coordinates.The lower panel measures the L error of this solution, testingentropy conservation (104), showing DISCO’s second-order con-vergence on this problem, including when the pulse is offset withrespect to the origin and mesh motion is utilized. for resolutions lower than 1024, convergence is fasterthan second order. At higher resolutions, convergenceis second order.This test is effectively 1D, as all motion is radial; it canbe run with arbitrarily low angular resolution. However,it can also be used as a multidimensional test problem, byoffsetting the origin of the pulse. Additional calculationswith these initial conditions were run with an origin offsetby ∆ y = 0 . Smooth Vortex
This test helps to demonstrate convergence, and the ef-fectiveness of the DISCO code at preservation of contactdiscontinuities. The grid has uniform radial resolutionfrom 0 < r <
5. Density is uniform with ρ = 1 .
0, andangular velocity is chosen as ω ( r ) = e − r (105)Pressure is chosen to balance centrifugal forces: ρ Ω r = ∂ r P . This results in the following pressure: P ( r ) = 1 − e − r . (106)4 Duffell Fig. 14.—
Passive scalar in the smooth vortex test at t = 10, demonstrating DISCO’s ability to maintain contact discontinuities to highprecision. The first two panels use N r = 64, with fixed mesh (left) and moving mesh (center). The fixed mesh diffuses out the passivescalar, whereas the moving mesh preserves the contact discontinuity precisely. A higher-resolution run ( N r = 256) is shown in the finalpanel for comparison. -7 -6 -5 -4 -2 L E rr o r Resolution Centered on origin
Fig. 15.—
Convergence of the smooth vortex test. Error isgiven by (107) and resolution indicates the number of radial zones.DISCO converges at second-order for this test.
The vortex is trans-sonic (maximum Mach number ofabout 0.53). In Figure 14, a passive scalar is included todemonstrate the code’s ability to maintain contact dis-continuities and prevent artificial diffusion. When meshmotion is turned off, the contact discontinuity is smearedout. With mesh motion turned on, the contact disconti-nuity is maintained precisely. The first two panels use alow resolution of 64 zones. A calculation at a resolutionof 256 radial zones is also included for comparison.Error is calculated at t = 10 using the density: L = R | ρ − . | dV R dV . (107)Figure 15 shows clear second-order convergence on thistest. Supersonic Keplerian Shear Flow
Fig. 16.—
Passive scalar in the Keplerian shear flow test, demon-strating DISCO’s ability to precisely preserve Keplerian orbitalflow, and to preserve contact discontinuities to high precision. Thepassive scalar is plotted after a single orbit at the outer boundary,corresponding to roughly 32 orbits at the inner boundary. Preserv-ing a Keplerian shear flow accurately is essential for any code whichis being used to study disks, in order that whatever phenomenonbeing studied is not swamped out by errors from this backgroundflow.
This is an important test, as most problems DISCOwas designed to solve have a Keplerian background flow.This stationary flow must be captured accurately if onewishes to study some subtle phenomenon which is a per-turbation to this flow.The setup for this problem is as follows. Zones arelogarithmically spaced in the range 0 . < r < .
0. Thedensity and pressure are uniform with ρ = 1 . P = 0 . v r = 0, ω ( r ) = r − / .Boundary conditions are fixed at these initial condi-tions. A point mass is inserted at r = 0 so that cen-trifugal and gravitational forces are balanced. For thisproblem, the gravitational potential was not smoothed,since the point mass is not on the grid, and the Keplerianflow is an exact solution. This results in a disk with aMach number of 7.7 at the outer boundary, and 24.5 atthe inner boundary.Because this test is axisymmetric (and therefore effec-ISCO Code 15 -5 -4 -3 -2 -2 L E rr o r Resolution
Fig. 17.—
Convergence of the Keplerian shear flow test. Error isgiven by (107) and resolution indicates the number of radial zones.DISCO achieves second-order convergence. tively 1D), it is not very important that the zones movewith the flow. However, for demonstration purposes apassive scalar has been added to the initial conditions: X = θ ( r cos( φ )) at t = 0. This passive scalar is plottedin Figure 16 at time t = 2 π , after the flow has had timeto shear it into a spiral.Error is computed identically to the vortex problem.This is computed at a time t = π , after a half-orbit at theouter boundary, and about 16 orbits at the inner bound-ary. Truncation error generates transient waves whichpropagate radially, bouncing between the two bound-aries. The L error is computed before these transientwaves have fully dissipated. Convergence is very close tosecond-order (Figure 17). For 512 radial zones, the L error is at the 10 − level.Because this problem is supersonic, most of the energyis kinetic, meaning that small relative errors in the en-ergy density can lead to large errors in the pressure. Forthis reason, it is important to evolve the modified energydensity, with Ω E ( r ) chosen to be Keplerian. Tests inwhich Ω E ( r ) was set to zero generated very large errorsnear the inner boundary. Cylindrical Kelvin-Helmholtz Instability
Flows unstable to Kelvin Helmholtz are traditionallytested on cartesian grids. However, Kelvin-Helmholtzinstability can occur in a rotational flow, as shown bythe following example. Radial zones are uniformly dis-tributed from 0 . < r < .
5. The background flow isgiven by a step function across r = 1:if r < ρ = 1 , ω = 2 , P = 4 + 2 r . (108)Otherwise, ρ = 2 , ω = 1 , P = 5 + r . (109)The perturbation is introduced in the radial velocity: Fig. 18.—
Density at t = 2 . v r = v cos (10 φ ) e − ( r − /σ (110)where v = 0 .
02 and σ = 0 .
1. Figure 18 plots the den-sity at time t = 2 .
0, showing that the instability is fullynonlinear at this time. The lower panels show the dif-ference when mesh motion is turned on and off. Sharperfeatures are present in the version in which the mesh ismoved, though care must be taken not to over-interprethow well these features are captured (Lecoanet et al.2016).
Viscosity
Cartesian Shear Flow with Viscosity | v | r cos( φ )Initial ConditionAnalytical SolutionNumerical Solution Fig. 19.—
Cartesian shear flow test, comparing DISCO’s solutionto the analytical Green’s function solution given by (112). Allviscous terms are tested here, as the flow is cartesian, and notaligned with the cylindrical grid.
In order to test DISCO’s implementation of viscosity, itis necessary to perform a noncircular viscosity test. Thisis so that every term in the viscous equations is used.Fortunately, this is as simple as setting up a cartesiantest problem and putting it on the cylindrical grid. Ofcourse, the cylindrical grid is not ideal for this test, butthe purpose is to make sure all of the terms (49) are im-plemented correctly, not to test accuracy or convergence.In cartesian coordinates, if we one designs a flow withuniform density and pressure and with ~v = v ( x )ˆ y , theNavier-Stokes equations reduce to:˙ v = νv ′′ , (111)where ν is viscosity. This has the well-known Green’sfunction solution: v y = v √ πνt exp { − ( x − x ) νt } (112)This solution is evolved using the following parameters: ρ = 1, P = 1, x = 1, v = 0 . ν = 0 . t = 1 .
0. 64 radialzones are uniformly distributed from 0 < r <
2. Figure19 plots the solution at this time as a function of the xcoordinate ( r cos( φ )) of each zone. Although the code isnot designed for problems so misaligned with the grid,the analytic solution is recovered.Additionally, the code was run with various terms inequations (49) set to zero, to check the importance ofeach term. The test has also been run with variouschoices of Ω E ( r ), to test that the solution is indepen-dent of this choice. All terms are necessary to capturethe solution to this level of accuracy. In other words, thistest ensures that the viscosity is implemented correctly. Supersonic Keplerian Spreading Test
A more complex viscosity test is attempted, relevantto accretion disks. A point mass is concentrated at r =0 and the domain extends from 0 . < r <
1. Orbitalvelocity is Keplerian, ω ( r ) = r − / , and pressure is setto a constant, P = 0 . ν = 10 − ensures an accretion flow given by v r = − ν/r (and Σ r Initial Condition1D SolutionDISCO Solution Fig. 20.—
Supersonic Keplerian spreading test. DISCO’s solu-tion is compared with the output of a 1D code which integratesthe 1D evolution equation for the surface density (114). the initial conditions assume this radial velocity). Thedensity is given by the following:Σ( r ) = 1 + 1 / √ r + e − r − . . (113)Note that “surface density” Σ is replacing “density” ρ ; this is a cosmetic change which reflects the fact thatthe code is integrating 1D and 2D (vertically integrated)disk equations.The first two terms are steady-state solutions to theevolution equations for Σ, but the final term is a densitybump which should be smeared out by viscosity. In thistest, a nearly isothermal equation of state is assumed( γ = 1 . r = 0 .
5. A very high Mach number is chosenbecause very cold disks are assumed when deriving the1D diffusion equation for the surface density:˙Σ = 3 r ( √ r ( √ r Σ ν ) ′ ) ′ (114)Similar to the previous test, analytic Green’s func-tion solutions exist to this equation. However, it is alsostraightforward to write a 1D code to integrate this PDEforward in time, to compare with DISCO’s solution. Forthe initial conditions stated above for the surface den-sity, a solution was obtained at t = 100 by integratingthis PDE.This is plotted in Figure 20, compared with DISCO’ssolution using 64 radial zones. The solutions do notmatch identically, but this is likely due to thin-disk as-sumptions made in deriving the 1D equation. Errors areat the percent-level. Disk-Planet Interactions
Low-Mass Planet
It is important to be able to capture planet-disk inter-actions in an idealized context. This test explores the lin-ear interaction between an Earth-mass planet and a su-personic Keplerian disk whose Mach number is M = 20at the orbital radius. The initial conditions are not meantto mimic a protoplanetary disk, rather they are meant toproduce an idealized environment in which it is easy toISCO Code 17 Fig. 21.—
Density for the Earth-mass planet in a Keplerian disk(after ten orbits, with 256 radial zones). The colormap is the sameas in Figure 12, but with density ranging between 0.97 and 1.03.Dashed curve is the analytical formula for the spiral wave given inequation (117). -1-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0 0.1 0 1 2 3 4 5 6 7 8 9 10 T o r que / T Time (orbits)-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 0.4 -3 -2 -1 0 1 2 3 T o r que D en s i t y [ q M ] (r-a)/h 64128256512 Fig. 22.—
Torque felt by the Earth-mass planet. The upperpanel plots torque density after ten orbits, showing that 128 radialzones is sufficient to capture this function. This is mirrored in thelower panel, which shows total torque as a function of time, atvarious resolutions. test the code. A small domain is used, with 0 . < r < . a = 1. Orbital veloc-ity is Keplerian ( ω ( r ) = r − / ), and surface density andpressure are uniform (Σ = 1, P = 0 . A ngu l a r M o m en t u m F l u x [ q M ] r/aRafikov (2002)12825651210242048 Fig. 23.—
Angular momentum flux carried by the planetarywake, for a planet with q = 10 − , corresponding to about threetimes Earth’s mass. The nonlinear propagation of the wave iscompared with the semi-analytic model of Rafikov (2002a). Thewave shocks a few scale heights away from the planet, outside ofwhich the angular momentum flux drops significantly. This weaknonlinear effect is much more difficult to capture than linear effectslike the torque density. P = Σ / M , (115)with boundary conditions fixed at the initial conditions.So far, this is just a supersonic Keplerian disk, similarto the test in section 3.1.4 (except with an isothermalequation of state). The only additional ingredient is apoint mass orbiting at a = 1.The planet mass is given by the Earth-to-Sun massratio: q = 3 × − = 0 . q NL (116)where q NL is defined according to the thermal massthreshold for nonlinearity ( q NL = M − = 1 . × − ).The planet’s gravity is given a smoothing length ǫ =0 . h = 0 . φ ( r ) = φ p + sign( r − a )(3 − p a/r − r/a ) M . (117)The spiral wave is established in the first orbit ( t = 2 π ).Figure 21 shows the density after 10 orbits, also plottingthe analytical prediction (117) for the spiral.Time-averaged torque density exerted by the planetis shown in Figure 22 for various resolutions, showingconvergence of the torque. The total torque is affectedby the nearby boundaries, but it reaches its convergedvalue with only 128 radial zones.Nonlinear propagation of the spiral wave is a muchmore subtle and difficult behavior to capture (Dong et al.2011). This can be measured by the angular momentumflux (AMF) emanating from the planet. This AMF hastwo components, a gravitational component due to theexcitation of the wave, and a “wave” component due tothe wave’s propagation and dissipation.Φ p = Φ grav + Φ wave (118)8 Duffell Fig. 24.—
Jupiter-mass planet in a Viscous Disk after 1000orbits, ν = 10 − . Because the planet mass is much larger than inprevious tests, the wave produces shocks strong enough to open adeep gap. Colormap is the same as in Figure 12, but with densityin logscale ranging from 10 − to 10 . Φ grav = R ∞ r dTdr dr r > a − R r dTdr dr r < a (119)Φ wave = Z rdφ Σ v r r ( ω − ω K ) , (120)where ω K is the background Keplerian orbital fre-quency. Goodman & Rafikov (2001) produced semi-analytical formulas for wave propagation and dissipation,which were later generalized to a global disk (Rafikov2002b). These scalings are plotted in Figure 23 com-pared with the measured planetary AMF for a planetwith q = 10 − (about three times Earth’s mass). Con-vergence is harder to establish for this test than for thelinear wave. 1024 zones (46 per scale height) are neededfor reasonably accurate measurement of the AMF. Viscous Disk with a Gap-Opening Planet
Larger planets produce stronger shocks, which dissi-pate the angular momentum in the spiral wave, deposit-ing torque in the disk, which causes an evacuation of anannulus in the vicinity of the planet’s orbit. This low-density annulus is known as a “gap”, and it is well-knownthat a Jupiter-mass planet can open a gap in a disk withmoderate viscosity.The scaling of gap depth with viscosity hasbeen determined in several numerical studies(Duffell & MacFadyen 2013; Fung et al. 2014;Kanagawa et al. 2015; Duffell 2015a):Σ gap Σ = 11 + f K ( q ) / π , (121)where f = 0 .
45 and K ( q ) ≡ q M /α. (122)This scaling will be measured using the same setupas the previous test, but with q = 10 − (Jupiter-to-sun Σ / Σ r 0.01 0.1 1 1 10 100 1000 Σ / Σ Time (orbits) ν = 3e-41e-43e-51e-53e-6 0.001 0.01 0.1 110 -6 -5 -4 -3 -2 Σ / Σ ViscosityNumerical Solution1D Model
Fig. 25.—
Gap depth for Jupiter in a viscous disk. The toppanel shows azimuthally averaged surface density as a function ofradius for the ν = 10 − case. The center panel shows gap depthas a function of time for various viscosities. The lower panel showsthe final gap depth (after 1000 orbits) as a function of viscosity,comparing with the analytical 1D model for the gap depth (121). mass ratio). Figure 24 shows the density at 10 orbitswith ν = 10 − . Figure 25 measures the gap depth asa function of time, and compares the steady-state gapdepth with the prediction from empirical calculationsand 1D models (121).Duffell & MacFadyen (2013) used this test as a meansto build a working definition of “numerical viscosity” inDISCO. The same test is performed with zero viscosity,and using a number of different numerical resolutions,and the gap depth was measured as a function of resolu-tion. This was compared to the gap depth as a functionof viscosity, to build a definition of “effective viscosity”of the numerical scheme. In that study, it was found that α num = 2 . × − (∆ r/h ) . (123) Magnetohydrodynamics
It will be important to test the accuracy and stabilityof the CT scheme described in section 2.10. Many workstesting numerical MHD studies look at the magnitude ofdiv B for given test problems, as a guide to code perfor-mance. Yet, the magnitude of div B can be small andstill affect the long-term behavior of the solution.The approach taken here to diagnosing div B perfor-mance is to measure the indirect impact of the diver-gence on the final solution. For example, one can advectISCO Code 19
Fig. 26.—
Magnetic field loop after advecting one orbit, using various methods. When constrained transport is turned off, B r is evolveddifferently from B φ , resulting in an asymmetric loop, and potentially resulting in unstable evolution. With CT turned on, and using theHLLD Riemann solver on the moving mesh, the loop can be maintained to machine precision. Colormap is the same as in Figure 12, butmagnetic energy is plotted, ranging from zero to 5 . × − . a field loop until ∇ · B errors produce inaccurate evolu-tion. Then one can demonstrate whether enforcing thediv B constraint resolves the issue. This approach leadsto a much greater confidence that attempts to eliminatemagnetic monopoles are accomplishing something. Orbital Advection of a Field Loop
The test which is typically most sensitive to div B er-rors is the advection of a magnetic field loop. A modi-fied version of this test is presented here to test DISCO’sCT algorithm in 2D. The initial setup (in the domain0 < r < .
5) is given as follows: ρ = 1 , ω = 1 , P = 0 .
01 + 0 . r . (124)The magnetic field is given by: B = (cid:26) B sin ( π ˜ r/R ) p r/R ˜ r < R r > R (125)where ˜ r is the distance from a point centered at r = 0 . φ = 0. R = 0 .
15 sets the radius of the loop. Thismagnetic field loop has magnetic tension and pressure,which can be balanced with the following adjustment tothe gas pressure: P → P − B ((˜ r/R )sin ( π ˜ r/R ) + (126)12 π ˜ r/R − π ˜ r/R ) + sin(4 π ˜ r/R )16 π )However, this adjustment is largely unnecessary, asthe magnetic field strength is very weak ( B = 10 − ).Stronger magnetic fields are possible to use, but theyare less susceptible to numerical disruption, as magnetictension and pressure provide restoring forces to the loop.Figure 26 shows the field loop after a single orbit ad-vected using N r = 128 radial zones, using 12 differentcombinations of schemes (with different options for theRiemann solver, mesh motion, and CT). HLLD shows asignificant improvement over the other schemes, reduc-ing the diffusion of the field loop. However, the loop ismore susceptible to div B instability when using this lessdiffusive Riemann solver. More diffusive methods miti-gate immediate catastrophe by diffusing out div B errors,but they do not eliminate the problem. Employing CTand HLLD, the loop is preserved to high precision, andwith the mesh motion the advection is solved to machineprecision (ignoring errors in the pressure balance, whichare small).0 Duffell Fig. 27.—
3D spinning field loop test after a single orbit, similarto the orbital advection test, but in three dimensions. A fixed meshsolution is compared with a moving mesh solution, showing thediffusion of the spinning loop that results from advection errors.Colormap is the same as in Figure 12, but magnetic energy isplotted, ranging from 0 to 5 . × − . Spinning 3D Loop
Testing the CT algorithm in 3D necessitates a 3D mag-netic field configuration, with nonzero components of B r , B φ and B z . The spinning 3D loop initially has field com-ponents B x and B z , and this loop rotates about the z-axisrigidly. Initial conditions are specified as: ρ = 1 . , ω = 1 . , P = 1 . . r . (127) B x = − B z/ ˜ r, B z = B x/ ˜ r (128)where ˜ r is a radial coordinate in the x-z plane:˜ r = p x + z (129) Fig. 28.—
2D MHD Flywheel test at t = 5. The upper panelhas constrained transport turned off, while the lower panel hasconstrained transport on. CT is necessary for numerical stabilityon this test; the grid-scale high-amplitude noise on the top panelis due to numerical instability. Colormap is the same as in Figure12, but magnetic energy is plotted, ranging from 0 to 4 . × − . B = ( B sin ( π ˜ r/R )cos( πy/R ) q rR ˜ r < R, | y | < . R B = 10 − (131)Figure 27 shows the spinning 3D loop using only 16radial zones (32 vertical zones) after a single orbit. Theupper panel uses a fixed mesh, while the lower panel usesa moving mesh. Moving the mesh allows for accuratepreservation of this loop. MHD Flywheel
The “MHD Flywheel” test is a new test presentedhere which illustrates the numerical instability associ-ated with div B errors. The setup involves a stationaryorbital MHD configuration, in which centrifugal force isISCO Code 21 -18 -16 -14 -12 -10 -8 -6 -4 -2
1 2 3 4 5 E rr o r ( C T O FF ) Time10 -18 -16 -14 -12 -10 -8 -6 -4 -2
1 2 3 4 5 E rr o r ( C T O N ) Time N R = 64128256512 Fig. 29.— L errors for the MHD flywheel test. Solid curvesindicate error in the density (equation 136), and dashed curves sig-nify error in the radial component of the magnetic field (equation137). The upper panel has uncontrolled divergence errors, whereasthe lower panel avoids this numerical instability by employing con-strained transport. balanced by magnetic tension, and magnetic pressure isbalanced by gas pressure: ρ = 1 , v r = 0 , ω = Ω e − r /R , (132) ~B = √ ρv φ ˆ φ, P = P − B . (133)This gives a stationary solution in 3D (though this testis only performed in 2D). Note that since velocity is par-allel to ~B , the total azimuthal flux of ~B is zero, whichmakes this test succeptible to div B violations if CT isnot employed. Constants are chosen so that P = 1 .
1( 12 e − ) ρ Ω R , (134)Ω = 1 . , R = 0 . . (135) Fig. 30.—
MHD explosion test at t = 0 .
2, using 256 radial zones.Colormap is the same as in Figure 12, but density ranges from 0 . .
7. The origin is not excised from the grid. DISCO is able tocapture all of the nontrivial shock structures with or without CT.
In 2D, this stationary solution is stable, so the goal isto preserve the initial conditions as precisely as possible.Figure 28 shows the solution at t = 5 using high resolu-tion ( N r = 512) and toggling CT on and off. While CTis turned off, div B errors cause a numerically unstablesolution. Implementing constrained transport eliminatesthese errors.This is shown quantitatively in Figure 29. Here, errorsin the solution are plotted. These are measured in boththe density and the radial component of the magneticfield: L ρ = R | ρ − . | dV R dV . (136) L B = R | B r | dV R dV . (137)Both errors are plotted in Figure 29, with CT turnedoff on the upper panel, and implemented in the lowerpanel. Without controlling div B, it is clear that er-rors in B r grow exponentially from machine round-off.Moreover, the growth rate appears to be proportionalto the resolution, so that improving the resolution givesgreater errors. Finally, when these errors grow from ma-chine round-off to be non-negligible contributions to thehydro evolution, they eventually affect the errors in den-sity (which would otherwise converge at second-order).By implementing constrained transport, the errors in B r remain very small. Cylindrical MHD Explosion
A commonly used setup which is useful for compar-isons with other codes, this test puts an explosion intoan initially uniform magnetic field, pushing the field linessideways and generating a nontrivial shock structure.The computational domain extends from 0 < r < . ρ = 1 .
0, and all velocities are initiallyzero. A magnetic field with strength B = 1 . Fig. 31.—
MHD rotor test at t = 0 .
15, using 256 radial zones.DISCO captures all of the details of the rotor, and is consistentwith other MHD codes on this test, whether or not CT is employed.Colormap is the same as in Figure 12, but gas pressure is plotted,ranging from 0 to 2. The origin is not excised from the grid. along the 45 ◦ diagonal (though the direction of the fieldis irrelevant, as the computational domain is cylindrical).The pressure is set to P = 0 . r < .
1, where P = 10. Figure 30 shows theexplosion at time t = 0 . N r = 256 radial zones(roughly equivalent to a cartesian box with 512 zonesacross). DISCO’s performance on this test can be com-pared e.g. with Figure 28 of Stone et al. (2008), or Figure4 of Pakmor et al. (2011).Note that the origin is kept on the grid, despite thecoordinate singularity at r = 0. In practice, the numer-ical solution to this test did not significantly depend onwhether CT was employed. MHD Rotor
This test is more challenging than the MHD explosion,and it tests the code’s propagation of strong torsionalAlfv´en waves. A rotational flow at the center winds upan initially uniform magnetic field: ω = ( v /R f ( r ) r < R v /rf ( r ) R < r < R r > R (138)where f ( r ) is a tapering function: f ( r ) = ( r < R ( R − r ) / ( R − R ) R < r < R r > R (139)Density, pressure, and radial velocity are given by ρ = 1 + 9 f ( r ) , P = 1 , v r = 0 (140)and the magnetic field is uniform in the x direction: ~B = (5 / √ π )ˆ x. (141)Another important point to note is that unlike the restof these test problems, which have either used γ = 5 / Fig. 32.—
3D MHD Explosion at time t = 0 .
1. Colormap is thesame as in Figure 12, but with density taking on values between 0and 2. This test uses 35 radial zones and 45 vertical zones. Thelower panel is the same test, but offset from the origin by x off = 0 . or an isothermal equation of state, this test uses a loweradiabatic index: γ = 1 . t = 0 .
15, which canbe compared with Figure 2 of Balsara & Spicer (1999) orFigure 25 of Stone et al. (2008).Again, the origin at r = 0 is resolved, and like theMHD explosion, the solution did not depend significantlyon whether CT was employed.
3D MHD Explosion
In order to test DISCO’s performance on a 3D MHDtest problem with nontrivial shocks, a 3D explosion testis performed. Initial conditions are given by ρ = 1 , ~B = 2ˆ z, (143) P = (cid:26) r < . . r > . t = 0 .
2. A second test was per-formed where the explosion is offset from the origin byan amount x = 0 .
1, in order to create a truly three-dimensional (non-axisymmetric) test; this is plotted inthe lower panel of Figure 32. This test is not performedin other studies, so there do not exist other examples toISCO Code 23
Fig. 33.—
Flock 3D MRI test after eight orbits, with N r × N z =512 ×
256 (x-z plane is shown). A coherent unstable mode emergesout of white noise initial perturbations (compare with Figure 5 ofFlock et al. 2010). Colormap is the same as in Figure 12, but withazimuthal magnetic field ranging from -0.02 to 0.02. -11 -10 -9 -8 -7 -6 -5 -4 -3 -2
2 4 6 8 10 R ad i a l M agne t i c E ne r g y Time (orbits) Nr x Nz = 32 x 1664 x 32128 x 64256 x 128512 x 2561024 x 512Flock et al (2010)
Fig. 34.—
Exponential growth of the MRI in the 3D Flock test.Dashed curve is the growth rate found by Flock et al. (2010) usingthe HLLD and Roe solvers using an upwind CT scheme. This rateis consistent with the analytical MRI growth rate. compare with, but the test demonstrates basic capturingof a 3D MHD shock structure, and that all geometricterms are correctly implemented, as the offset test hasan equivalent solution to the centered test. Additionally,this test shows that DISCO is capable of evolving shockspropagating along the coordinate axis, and shocks collid-ing with the coordinate singularity at r = 0, especiallyin the offset case, where the pressure jump is initiatedacross r = 0 in the equatorial plane. Magnetorotational Instability: Linear Growth
One of the most important astrophysical applicationsof numerical MHD is the study of the magnetorotationalinstability (MRI Balbus & Hawley 1991; Hawley et al.1995; Stone et al. 1996). MRI is thought to be the sourceof angular momentum transport (and hence accretion) inmost astrophysical disks. Growth of MRI using DISCOwill be studied using two test problems. First, lineargrowth is explored, by implementing the 3D MRI test
Fig. 35.—
Nonlinear 3D MRI test after 25 orbits. Logarithmof magnetic energy is plotted. This test demonstrates DISCO’sability to capture MRI, not just in the linear growth phase, but inthe nonlinear, fully turbulent regime. of Flock et al. (2010). The set-up was described veryclearly in that study, so the initial conditions will not berepeated here.The velocity field is seeded with white noise, and themagnetic field is chosen so that n = 4 corresponds tothe fastest-growing MRI mode. Figure 33 shows the az-imuthal component of the magnetic field after eight or-bits, to be compared with Figure 5 of Flock et al. (2010).Coherent magnetic fields have grown from white noiseinitial perturbations. Growth of the MRI is shown inFigure 34, which plots the “Radial magnetic energy”,showing the exponential growth of the instability, com-pared with the growth rate found by Flock et al. (2010)(see their Figure 3). At resolutions above N z = 32, or8 zones per MRI wavelength, the correct growth rate isrecovered. Magnetorotational Instability: Nonlinear Saturation
Capturing the correct linear growth rate of MRI isan important benchmark for any 3D MHD code, but itshould be recognized that linear growth does not con-stitute a complete MRI test. The true test comes inthe nonlinear, fully turbulent phase, when the statisti-cal properties of the turbulence could be affected by thenumerical scheme.The nonlinear MRI test of Sorathia et al. (2012) is im-plemented here, to determine DISCO’s performance ona 3D nonlinear turbulent magneto-rotational flow. Theinitial conditions used are given by their “zero net flux”case. Figure 35 shows a snapshot in time of the magneticenergy, to be compared with Figure 1 of Sorathia et al.(2012). Figure 36 plots the quantities α M and β , whichare given by4 Duffell α M Zones per scale height = 8163264 ( π /4 Wedge)0.0000.0200.0400.0600.080 β - θ B ( deg r ee s ) Time (Orbits)
Fig. 36.—
Growth and saturation in the 3D nonlinear MRI test.The various resolutions shown covered two scale heights, so that N r × N z = 64 ×
16, 128 ×
32, 256 ×
64, and 512 × π in azimuth. Thecyan curve is a calculation on a wedge spanning π/ ◦ , consistent with the anglefound in the present study (bottom panel). α M = − h B r B φ i / h P i , (145) β = h P i / (cid:28) B (cid:29) . (146)where the average is taken over the entire domain. Thesecan be compared with Sorathia et al. (2012), Figure 5.Finally, the magnetic “Tilt angle” is measured. This isdefined as θ B = sin − ( αβ ) / . (147)Sorathia et al. (2012) showed that this tilt angle is ro-bustly 12 ◦ (see their Figure 11 and their Table 3). Inthe bottom panel of Figure 36, this tilt angle is plotted,showing consistency with this value. Efficiency and Scaling
Finally, a few tests are performed to determineDISCO’s performance and scaling to large numbers ofprocessors. The Keplerian shear flow test (section 3.1.4)is performed in both 2D and 3D, on various num-bers of processors, and with variable resolution. Both“strong scaling” (with fixed resolution) and “weak scal-ing” (with fixed number of zones per process) are mea-sured in Figure 37. Calculations wer performed with 2.8GHz Intel Xeon E5-2680v2 processors, using the 20-core“Ivy-Bridge” nodes on NASA’s Pleiades supercomputer.
1 10 100 1000 10000 Z one upda t e s pe r C P U s e c ond Number of CPUsWeak Scaling, Zones / CPU = 1.4 x 10 Strong Scaling, N R = 4096 Fig. 37.—
Efficiency and scaling of DISCO for the Keplerianshear flow test (section 3.1.4). When run in parallel, DISCO typi-cally achieves about 5 × zone updates per CPU second. Bothstrong and weak scaling preserve this rate up to thousands of CPUs. DISCO’s peak performance is 9 × zone updates perCPU second. In the range between a few nodes and up tothousands of CPUs, typical performance is 5 × zoneupdates per CPU second. A latency is reached whenthere are as few as two radial zones per CPU in 2D (ortwo radial and two vertical zones per CPU in 3D). DISCUSSION
A moving-mesh technique is presented for numeri-cally integrating both hydro and MHD equations in3D, with specific application to astrophysical disks,using a new constrained transport technique. TheDISCO code has been made publicly available at https://github.com/duffell/Disco under the GNUgeneral public license.The numerical scheme has been detailed in this work,including a description of the novel constrained transporttechnique. The versatility of the scheme is shown by de-scribing the implementation of many different systems ofequations, including the “bare” Euler equations, viscoushydrodynamics, magnetohydrodynamics, and includingadditional terms due to the gravitational influence of or-biting bodies.Many code tests have been performed, demonstratingDISCO’s ability to integrate all of these different systemsof equations. DISCO has no trouble accurately captur-ing shocks, whether they are aligned or misaligned withthe orbital motion. DISCO excels at advecting contactdiscontinuities and MHD discontinuities with the orbitalflow. DISCO also converges at second-order for smoothflows, and can evolve flows with very high Mach numberwithout accumulating significant errors.DISCO is ideal for studies of disk-satellite interactions,especially problems which involve a near-equal-mass bi-nary, where both point masses must be resolved on thegrid. It is also an accurate code for calculating MRI, asdemonstrated by code tests of both linear stability andnonlinear turbulence.The constrained transport technique for maintainingzero divergence errors while updating the magnetic fieldsis unique and applicable to nontrivial and dynamic meshtopologies. Stable, accurate evolution is possible withthis CT method and, as demonstrated in several of theISCO Code 25code tests, CT is necessary to prevent inaccurate or un-stable evolution.Certainty in the robustness of the CT method is as-sured by DISCO’s accurate calculation of the magnetoro-tational instability, capturing both the linear growth rateand the nonlinear turbulent phase in quantitative com-parisons with previous numerical studies.In the future, it will also be possible to use DISCO tointegrate the equations of general relativistic magneto-hydrodynamics (GRMHD), for application to black holeaccretion disks. Self-gravity is another future improve-ment which is possible. Both of these improvements willbe presented in a future study.Resources supporting this work were provided by theNASA High-End Computing (HEC) Program throughthe NASA Advanced Supercomputing (NAS) Divisionat Ames Research Center. Some early exploratory cal-culations employed the Savio computational cluster re-source provided by the Berkeley Research Computingprogram at the University of California, Berkeley (sup-ported by the UC Berkeley Chancellor, Vice Chancellorof Research, and Office of the CIO). I am grateful to An-drew MacFadyen, Daniel D’Orazio, Brian Farris, JeffreyFung, Andrei Gruzinov, Zoltan Haiman, Frederic Mas-set, Philip Mocz, Diego Mu˜noz, Eliot Quataert, StephanRosswog, Geoff Ryan, Jim Stone, and Jonathan Zrake forhelpful comments and discussions. I would also like tothank the anonymous referee for the very helpful review.
APPENDIX
VISCOSITY
In section 2.9, viscous terms were added to Euler’sequations. Here these terms are derived. The Navier-Stokes equation is given by the following: ∂ t ( ρ~v ) + ∇ · ( ρ~v~v ) + ~ ∇ P = ∇ · σ, (A1)Where σ is the viscous stress tensor: σ ↔ = ρν ( ~ ∇ ~v + ~ ∇ ~v + η ∇ · v I ↔ ) (A2)The quantity η is related to the ratio of bulk viscosityto shear viscosity. For pure shear flow, the value of η should have no impact on the solution. Initially, the valueof η = 0 is chosen, and will be re-introduced near the endof the derivation. Additionally, for ease of presentation,the prefactor ρν will be ignored (set to unity) until thelast few steps. Finally, this derivation will be carried outin 2D ( r, φ ) for simplicity, but the additional terms inthe vertical dimension are straightforward to calculate.The derivatives in the formula for σ are covariantderivatives, so in index notation the expression evaluatesas: σ ij = ∇ i v j + ∇ j v i = ∂ i v j + ∂ j v i − kij v k , (A3)where Γ is the Christoffel symbol for flat space in cylin-drical coordinates: Γ rφφ = − r (A4)Γ φφr = Γ φrφ = 1 /r, (A5) and all other components evaluate to zero. Notation willbe used which removes the confusion with indices, in-troducing the nonindexed v for radial velocity and ω forangular velocity: v r = v r ≡ v, (A6) v φ = ω, v φ = r ω (A7)Also, for brevity, partial derivatives in radius are re-placed with primes: σ φφ = 2 ∂ φ ( r ω ) + 2 rv (A8) σ φr = ∂ φ v + r ω ′ (A9) σ rr = 2 v ′ (A10)Note that these expressions for the viscous stress agreewith others in the literature, except that these formulasare expressed in a coordinate basis, whereas often anorthonormal basis is used: σ ˆ φ ˆ φ = 2 r ( ∂ φ v ˆ φ + v ˆ r ) (A11) σ ˆ φ ˆ r = ∂ r v ˆ φ + 1 r ( ∂ φ v ˆ r − v ˆ φ ) (A12) σ ˆ r ˆ r = 2 ∂ r v ˆ r (A13)This result agrees, for example, with Edgar (2006),equations (10), (11) and (13), except of course thatEdgar’s result does not set η = 0.However, this result does not agree at face value withstandard expressions in for viscous terms in cylindricalcoordinates (e.g. Landau & Lifshitz 1966, page 51 of thefirst edition). It is not enough to evaluate the differentcomponents of σ ↔ in cylindrical coordinates, since thereare additional geometric terms which come from takingits divergence. If σ ↔ were a vector, its components wouldbe sufficient, since one can use the divergence theoremto evaluate it. However, since σ is a tensor, additionalterms appear. The appropriate goal, then, is to write σ in the following form:( ∇ · σ ) j = 1 √ g ∂ r ( √ g ˜ σ rj ) + 1 √ g ∇ φ ( √ g ˜ σ φj ) + S j (A14)This expression looks like the divergence of a vectorplus a source term, and therefore ˜ σ ij can be manipu-lated as a vector (treating the raised index as a vectorindex). It is important to point out that the presence ofa source term here is inescapable, due to the fact thatthe geometry does not have a symmetry in the radialdirection.It is straightforward to express the divergence of a ten-sor in this form, using simple Riemannian geometry:( ∇ · σ ) j = 1 √ g ∂ i ( √ gσ ij ) − ∂ j g kl σ kl (A15)Looking at the radial component:( ∇ · σ ) r = 1 r ∂ r ( rσ rr ) + 1 r ∂ φ ( rσ φr ) − rσ φφ (A16)6 Duffell= 1 √ g ∂ r ( √ g v ′ )+ 1 √ g ∇ φ ( √ g ( ∇ φ v + rω ′ − ω )) (A17) − v/r Similarly for the azimuthal component:( ∇ · σ ) φ = 1 r ∂ r ( rσ rφ ) + 1 r ∂ φ ( rσ φφ ) (A18)= 1 √ g ∂ r ( √ g ( r ω ′ + ∂ φ v ))+ 1 √ g ∇ φ ( √ g (2 r ∇ φ ω + 2 v )) (A19)There are several things to note about these ex-pressions. First, they also don’t agree with stan-dard formulas for viscosity in cylindrical coordinates(Landau & Lifshitz 1966). In particular, there are mixedderivative terms. This can be resolved by choosing asuitable value of η . Recall that η is the coefficient of thedivergence of velocity, which will be abbreviated as θ : θ = ∇ · v = v ′ + v/r + ∂ φ ω (A20)The η term modifies the viscous tress as: σ ij → σ ij + ηθg ij (A21)Component-wise: σ rr → σ rr + ηθ (A22) σ φφ → σ φφ + r ηθ (A23)This has the following impact on the divergence: δ ( ∇ · σ ) r = 1 r ∂ r ( rηθ ) − r ηθ = ∂ r ( ηθ ) (A24)= η (cid:18) r ∂ r ( rv ′ ) − vr + 1 r ∂ φ ( rω ′ ) (cid:19) . (A25)Similarly, for the φ component: δ ( ∇ · σ ) φ = η (cid:18) r ∂ r ( r∂ φ v ) + ∇ φ ( r ∇ φ ω ) (cid:19) (A26)There is a natural choice for η , η = − , (A27)which simplifies all of these formulas. This eliminatesall of the mixed derivative terms, and restores agreementwith standard formulas for viscosity in cylindrical coor-dinates: ( ∇ · σ ) r = 1 √ g ∂ r ( √ gv ′ )+ 1 √ g ∇ φ ( √ g ( ∇ φ v − ω )) (A28) − v/r ( ∇ · σ ) φ = 1 √ g ∂ r ( √ gr ω ′ )+ 1 √ g ∇ φ ( √ g ( r ∇ φ ω + 2 v )) (A29)This finally results in the following ˜ σ :˜ σ rr = v ′ (A30)˜ σ rφ = ∇ φ v − ω (A31)˜ σ φr = r ω ′ (A32)˜ σ φφ = r ∇ φ ω + 2 v (A33)And a source term for the momentum in the radialdirection: S r = − v/r (A34)As expected, there is no source term for angular mo-mentum. After re-introducing the coefficient ρν , equa-tions (A30-A34) can be re-expressed as the terms givenin (49) and (50).As mentioned in Section 2.9, introducing Ω E ( r ) pro-duces a source term for the energy, given by S Energy visc = σ rφ r Ω ′ E ( r ) . (A35)using the derived formula for σ , this can be expressed as S Energy visc = ρν ( ∇ φ v + rω ′ ) r Ω ′ E ( r ) . (A36) REFERENCESAbbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, PhysicalReview Letters, 116, 061102Abramowicz, M. A., & Fragile, P. C. 2013, Living Reviews inRelativity, 16, arXiv:1104.5499Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza,L. 2014, Protostars and Planets VI, 475Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Balsara, D. S., & Spicer, D. S. 1999, Journal of ComputationalPhysics, 149, 270Ben´ıtez-Llambay, P., & Masset, F. 2016, ArXiv e-prints,arXiv:1602.02359 Brackbill, J. U., & Barnes, D. C. 1980, Journal of ComputationalPhysics, 35, 426Brio, M., & Wu, C. C. 1988, Journal of Computational Physics,75, 400Dai, W., & Woodward, P. R. 1998, ApJ, 494, 317DeBuhr, J., Zhang, B., Anderson, M., Neilsen, D., & Hirschmann,E. W. 2015, ArXiv e-prints, arXiv:1512.00386Dedner, A., Kemm, F., Kr¨oner, D., et al. 2002, Journal ofComputational Physics, 175, 645Del Zanna, L., Bucciantini, N., & Londrillo, P. 2003, A&A, 400,397