VADER: A Flexible, Robust, Open-Source Code for Simulating Viscous Thin Accretion Disks
VVADER : A Flexible, Robust, Open-Source Code forSimulating Viscous Thin Accretion Disks
Mark R. Krumholz and John C. Forbes
Astronomy Department, University of California, Santa Cruz, 95064, USA;[email protected] and [email protected]
Abstract
The evolution of thin axisymmetric viscous accretion disks is a classic prob-lem in astrophysics. While models based on this simplified geometry provideonly approximations to the true processes of instability-driven mass and an-gular momentum transport, their simplicity makes them invaluable tools forboth semi-analytic modeling and simulations of long-term evolution wheretwo- or three-dimensional calculations are too computationally costly. De-spite the utility of these models, the only publicly-available frameworks forsimulating them are rather specialized and non-general. Here we describe ahighly flexible, general numerical method for simulating viscous thin diskswith arbitrary rotation curves, viscosities, boundary conditions, grid spac-ings, equations of state, and rates of gain or loss of mass (e.g., through winds)and energy (e.g., through radiation). Our method is based on a conservative,finite-volume, second-order accurate discretization of the equations, which wesolve using an unconditionally-stable implicit scheme. We implement Ander-son acceleration to speed convergence of the scheme, and show that thisleads to factor of ∼ VADER ), which is freely available for download from https://bitbucket.org/krumholz/vader/ under the terms of the GNUGeneral Public License.
Keywords: accretion, accretion disks — applied computing ∼ astronomy —mathematics of computing ∼ partial differential equations — mathematics ofcomputing ∼ nonlinear equations — methods: numerical Preprint submitted to Astronomy & Computing February 25, 2015 a r X i v : . [ a s t r o - ph . I M ] F e b . Introduction Accretion disks are ubiquitous in astrophysics, in fields ranging from starand planet formation to high energy astrophysics to galaxies, and an enor-mous amount of effort has been invested in modeling them (e.g. Pringle,1981). One approach to constructing such models is to conduct full two- orthree-dimensional simulations, and this method offers the highest fidelity tothe actual physical processes taking place in disks. However, such simula-tions are impractically computationally expensive for phenomena that takeplace over very large numbers of orbital timescales, or for disks where thecharacteristic scales that must be resolved for the simulation to convergeare vastly smaller than the disk radial extent. For example, a long-term two-dimensional simulation of a protoplanetary disk might cover several thousandorbits, but the timescale over which planets form is millions of orbits. Quasi-periodic oscillations from disks around black holes, neutron stars, and whitedwarfs can take place over similarly large numbers of orbits. In galaxies,the number of orbits is relatively modest, but the characteristic size scale ofgravitational instability for the coldest phase of the interstellar medium is ∼ times smaller than the disk radial extent. None of these problems areamenable to solution by two- or three-dimensional simulations, at least notwithout extensive use of sub-grid models to ease the resolution requirements.In such cases, one-dimensional simulations in which the disk is treated asvertically thin and axisymmetric are a standard modeling tool. The generalapproach in such simulations is to approximate the turbulence responsiblefor transporting mass and angular momentum through the disk as a vis-cosity, and to develop an analytic or semi-analytic model for this transportmechanism. Cast in this form, the evolution of a disk is described by a pairof one-dimensional parabolic partial differential equations for the transportof mass and energy; the form of these equations is analogous to a diffusionequation in cylindrical coordinates. Depending on the nature of the problem,these equations may have source or sink terms, may have a wide range ofboundary conditions, and may have multiple sources of non-linearity.Thus far in the astrophysical community most viscous disk evolutioncodes have been single-purpose, intended for particular physical regimes andmodeling physical processes relevant to that regime. Thus for example thereare codes intended for protoplanetary disks that include models for accretiononto the disk during ongoing collapse (e.g. Hueso and Guillot, 2005; Visserand Dullemond, 2010; Lyra et al., 2010; Horn et al., 2012; Benz et al., 2014),2alaxy disk codes containing prescriptions for star formation (e.g. Forbeset al., 2012, 2014), and codes for simulating accretion onto compact objectsthat include models for magnetically-dominated coronae and have equationsof state that include the radiation pressure-dominated regime (e.g. Liu et al.,2002; Mayer and Pringle, 2007; Cambier and Smith, 2013). While these codesare specialized to their particular problems, they are often solving very simi-lar systems of equations, and thus there is a great deal of replication of effortin every community developing its own code.This is particularly true because most of the codes are not open source,and for the most part the authors have not published detailed descriptionsof their methodologies, forcing others to invent or re-discover their own. Thesole exceptions of which we are aware are the GIDGET code for simulatinggalaxy disk evolution (Forbes et al., 2012) and the α -disk code publishedby Lyra et al. (2010) and Horn et al. (2012), based on the PENCIL code(Brandenburg and Dobler, 2002). Neither
GIDGET nor Lyra’s code are suitedfor general use. For example, neither allows a wide range of viscosities,equations of state, and rotation curves. Ironically, this situation is in sharpcontrast to the situation for two-dimensional disk simulations, where thereare a number of open source codes that include viscosity and various otherphysical processes. These include
ZEUS-2D (Stone and Norman, 1992a,b;Stone et al., 1992),
VHD (McKinney and Gammie, 2002, 2013), and
PLUTO (Mignone et al., 2007). However, while it is possible to run all these codesin either a one-dimensional or pseudo-one-dimensional mode, they are allbased on explicit schemes, which limits their ability to simulate very longtime scales.The goal of this paper is to introduce a very general method for computingthe time evolution of viscous, thin, axisymmetric disks in one dimension,using a method suitable for simulating disks over many viscous evolutiontimes at modest computational cost. We embed this method in a code calledViscous Accretion Disk Evolution Resource (
VADER ), which we have releasedunder the GNU General Public License.
VADER is available for download from https://bitbucket.org/krumholz/vader/ . The code is highly flexible andmodular, and allows users to specify arbitrary rotation curves, equations ofstate, prescriptions for the viscosity, grid geometries, boundary conditions,and source terms for both mass and energy. The equations are written inconservation form, and the resulting algorithm conserves mass, momentum,angular momentum, and energy to machine precision, as is highly desirablefor simulations of very long term evolution. We employ an implicit numerical3ethod that is unconditionally stable, allows very large time steps, and is fastthanks to modern convergence acceleration techniques.
VADER is descendedfrom
GIDGET (Forbes et al., 2012, 2014) in a very general sense, but it isdesigned to be much more flexible and modular, while omitting many of thefeatures (e.g., cosmological accretion and the dynamics of collisionless stars)that are specific to the problem of galaxy formation. It is implemented in Cand Python. The present version is written for single processors, but we planto develop a threaded version in the future once an open source, threadedtridiagonal matrix solver becomes available.The plan for the remainder of this paper is as follows. In Section 2we introduce the underlying equations that
VADER solves, and describe ouralgorithm for solving them. In Section 3 we present a number of tests ofthe code’s accuracy and convergence characteristics. Section 4 discusses theefficiency and performance of the algorithm. Finally we summarize in Section5.
2. Equations and Simulation Algorithm
The physical system that
VADER models is a thin, axisymmetric disk ofmaterial in a time-steady gravitational potential. We consider such a diskcentered at the origin and lying in the z = 0 plane of a cylindrical ( r, φ, z )coordinate system. The equations of continuity and total energy conservationfor such a system, written in conservation form, are (e.g., equations 1 andA13 of Krumholz and Burkert 2010) ∂∂t Σ + 1 r ∂∂r ( rv r Σ) = ˙Σ src (1) ∂∂t E + 1 r ∂∂r [ rv r ( E + P )] − r ∂∂r (cid:18) r v φ T πr (cid:19) = ˙ E src . (2)Here Σ is the mass surface density in the disk, E = Σ (cid:18) v φ ψ (cid:19) + E int ≡ Σ ψ eff + E int (3)is the total energy per unit area, v φ is the rotation speed as a function ofradius, ψ is the gravitational potential (which is related to v φ by ∂ψ/∂r = v φ /r ), ψ eff = ψ + v φ / int is the internal energy per unit area, P is the vertically-integrated pressure( (cid:82) ∞−∞ p dz , where p is the pressure), v r is the radial velocity, and T is thetorque applied by a ring of material at radius r to the adjacent ring at r + dr . The source terms ˙Σ src and ˙ E src represent changes in the local massand energy per unit area due to vertical transport of mass (e.g., accretionfrom above, mass loss due to winds, or transformation of gas into collisionlessstars) or energy (e.g., radiative heating or cooling). VADER allows very generalequations of state; the vertically-integrated pressure P may be an arbitraryfunction of r , Σ, and E int (but not an explicit function of time).The torque and radial velocity are related via angular momentum con-servation, which implies v r = 12 πr Σ v φ (1 + β ) ∂∂r T , (4)where β = ∂ ln v φ /∂ ln r is the index of the rotation curve at radius r . To pro-ceed further we require a closure relation for the torque T . For the purposesof this calculation we adopt the Shakura and Sunyaev (1973) parameteriza-tion of this relation, slightly modified as proposed by Shu (1992). In thisparameterization, the viscosity is described by a dimensionless parameter α such that T = − πr α (1 − β ) P. (5)Note that inclusion of the 1 − β term is the modification proposed by Shu(1992), and simply serves to ensure that the torque remains proportional tothe local rate of shear in a disk with constant α but non-constant β . Withthis definition of α , the kinematic viscosity is ν = α (cid:18) rv φ (cid:19) (cid:18) P Σ (cid:19) . (6)The dimensionless viscosity α (as well as the source terms ˙Σ src and ˙ E src ) canvary arbitrarily with position, time, Σ, and E .Since these equations are derived in Krumholz and Burkert (2010), we willnot re-derive them here, but we will pause to comment on the assumptionsthat underlie them, and the potential limitations those assumptions imply.The system of equations is appropriate for a slowly-evolving thin disk withnegligible radial pressure support. Specifically, we assume that (1) the scaleheight H (cid:28) r (thinness), (2) the radial velocity v r obeys v r (cid:28) v φ and5 v r (cid:28) E int (slow evolution), and (3) Σ | v φ /r − ∂ψ/∂r | (cid:28) E int (negligibleradial pressure support). VADER is not appropriate for disks that do notsatisfy these assumptions.
To discretize the equation in space, consider a grid of N cells with edges lo-cated at positions r − / , r / , r / , . . . r N − / and centers at positions r , r , r , . . . r N − .Most often the grid will be uniformly-spaced in either the logarithm of r , inwhich case r i = √ r i − / r i +1 / and ∆ ln r i +1 / = ln( r i +1 /r i ) is constant, or in r itself, so that r i = ( r i − / + r i +1 / ) / r i +1 / = r i +1 − r i is constant.However, VADER allows arbitrary placement of the cell edges, so long as thesequence r i +1 / is strictly increasing with i .Let A i = π ( r i +1 / − r i − / ) be the area of cell i . We will similarly denotethe rotation curve and its logarithmic derivative evaluated at cell centersand edges by v φ,i , β i , v φ,i +1 / , and β i +1 / . Integrating equations (1) and (2)over the area of cell i , and making use of the divergence theorem to evaluateterms involving the operator (1 /r )( ∂/∂r )( r · ) (which is simply the radialcomponent of the divergence operator written out in cylindrical coordinates)gives ∂∂t Σ i + F M,i +1 / − F M,i − / A i = ˙Σ src ,i (7) ∂∂t E i + F E,i +1 / − F E,i − / A i + F T,i +1 / − F T,i − / A i = ˙ E src ,i , (8)where we have defined the surface density Σ i averaged over cell i byΣ i = 1 A i (cid:90) A i Σ dA, (9)and similarly for E i , ˙Σ src ,i , and ˙ E src ,i , and we have defined the fluxes at celledges by F M,i +1 / = (2 πrv r Σ) r = r i +1 / (10) F E,i +1 / = [2 πrv r ( E + P )] r = r i +1 / (11) F T,i +1 / = [2 πrv φ α (1 − β ) P ] r = r i +1 / , (12)and similarly for i − /
2. The first two terms above represent the advectivefluxes of mass and enthalpy (including gravitational potential energy in the6nthalpy), respectively, while the third term is the energy flux associatedwith work done by one ring on its neighbors through the viscous torque.Note that the fluxes as defined here are total fluxes with units of mass orenergy per time, rather than flux densities with units of mass or energy pertime per area.For the purposes of numerical computation it is convenient to replace theenergy equation with one for the pressure. We let E i = E int ( r, Σ i , P i ) + Σ i ψ eff ,i , (13)where E int ( r, Σ , P ) is the function giving the relationship between internalenergy per unit area, surface density, and vertically-integrated pressure. Notethat, for an ideal gas, E int is a function of P alone, and if the disk is vertically-isothermal then it is given by E int = P/ ( γ − γ is the ratio of specificheats. However, we retain the general case where γ can be an explicit functionof Σ, P , and r (but not time) to allow for more complex equations of state.We do not, however, allow γ = 1 exactly, because in the truly isothermal casethe energy equation vanishes and the character of the system to be solvedchanges fundamentally. We show below that one can approximate isothermalbehavior simply by setting γ = 1 + (cid:15) , where (cid:15) is a very small parameter, andthat with this approach VADER has no difficulty recovering analytic solutionsthat apply to truly isothermal disks.Substituting this form for E i in equation (8), applying the chain rule, andusing equation (7) to eliminate ∂ Σ i /∂t gives ∂∂t P i + F P,i + − F P,i − A i = ( γ i −
1) ˙ E int , src ,i . (14)Here γ ≡ ∂P∂E int (cid:12)(cid:12)(cid:12)(cid:12) r, Σ , (15)and the source term ˙ E int , src ≡ ˙ E src − (cid:18) ψ eff + δ P Σ (cid:19) ˙Σ src , (16)where δ ≡ γ − ∂ ln P∂ ln Σ (cid:12)(cid:12)(cid:12)(cid:12) r,E int . (17)7or a vertically-isothermal ideal gas, δ = 0 because P does not change asΣ is varied at fixed E int . From the definition of ˙ E int , src we can see that thisterm represents the time rate of change of the internal energy due to externalforcing (e.g., radiative losses) evaluated at fixed Σ. For example, if a diskis both cooling radiatively and losing mass to a wind, then ˙ E int , src includesthe rate of change of the internal energy due to radiation alone, but not anychange in the internal energy due to mass loss from the wind. Finally, wehave defined the left and right pressure fluxes in cell i by F P,i − = ( γ i − (cid:20) F E,i − / + F T,i − / − (cid:18) ψ eff ,i + δ i P i Σ i (cid:19) F M,i − / (cid:21) (18) F P,i + = ( γ i − (cid:20) F E,i +1 / + F T,i +1 / − (cid:18) ψ eff ,i + δ i P i Σ i (cid:19) F M,i +1 / (cid:21) . (19)Note that in general F P,i + (cid:54) = F P, ( i +1) − . An important point to emphasize isthat equation (14) is derived from the already-discretized energy equation(8). Thus, if γ and δ are known analytically, which is the case for anysimple equation of state, then a numerical implementation of equation (14)will conserve energy to machine precision. If γ and δ must be obtainednumerically, then the accuracy of conservation would depend on the accuracywith which they are determined. In fact, as we discuss below, when γ and δ are non-constant we use a scheme that is not explicitly conservative inany event, though in practical tests we find that the accuracy of energyconservation is within ∼ v r , which in turn depends on the partial derivative ofthe torque and thus of the pressure. We approximate this using centereddifferences evaluated at the cell edges. Even if the actual distribution ofcell positions is uniform in neither linear or logarithmic radius, we requirethat a grid be classified as log-like or linear-like. For a linear grid we usethe standard second-order accurate centered difference, while for logarithmicgrids we rewrite the derivative ∂ T /∂r as (1 /r ) ∂ T /∂ ln r and use a centereddifference to evaluate the derivative with respect to ln r . This guarantees thatthe derivatives are second-order accurate in either r or ln r for uniformly-spaced linear or logarithmic grids; for non-uniform grids the derivatives arefirst-order accurate. This choice gives a mass flux F M,i +1 / = − g i +1 / (cid:2) α i +1 (1 − β i +1 ) r i +1 P i +1 − α i (1 − β i ) r i P i (cid:3) , (20)8here for convenience we have defined the factor g i +1 / = 2 πv φ,i +1 / (1 + β i +1 / ) (cid:26) / ( r i +1 / ∆ ln r i +1 / ) , (logarithmic grid)1 / ∆ r i +1 / , (linear grid) (21)for logarithmic and linear grids, respectively. Using the same strategy forthe torque flux yields F T,i +1 / = πr i +1 / v φ,i +1 / (1 − β i +1 / ) ( α i +1 P i +1 + α i P i ) . (22)Finally, the enthalpy flux is F E,i +1 / = h i +1 / F M,i +1 / , (23)where h i +1 / = [( E + P ) / Σ] i +1 / is our estimate for the specific enthalpy(including the gravitational potential energy) at the cell edge. To estimate h i +1 / , we divide the enthalpy into an internal part and a gravitational plusorbital part, i.e., we set h i +1 / = (cid:18) E int + P Σ (cid:19) i +1 / + ψ eff ,i +1 / (24) ≡ h int ,i +1 / + ψ eff ,i +1 / . (25)The gravitational plus orbital part ψ eff ,i +1 / is known exactly from the spec-ification of the rotation curve, so no approximation is necessary for it. Toestimate h int ,i +1 / , VADER uses a first-order upwind scheme (Fletcher, 1991)to maintain stability: h int ,i +1 / = F M,i +1 / + h int ,i +1 / ,L + F M,i +1 / − h int ,i +1 / ,R (26)where F M,i +1 / + = max( F M,i +1 / ,
0) (27) F M,i +1 / − = min( F M,i +1 / , . (28)The left and right internal enthalpies h int ,i +1 / ,L and h int ,i +1 / ,R represent thevalues on the left and right sides of the interface. VADER can set these valuesusing either piecewise constant, slope-limited piecewise linear, or piecewiseparabolic extrapolation, yielding first-order accurate, second-order accurate,and third-order accurate approximations, respectively. If h int ,i = [( E int +9 ) / Σ] i is the internal enthalpy evaluated at the cell center, then piecewiseconstant interpolation gives h int ,i +1 / ,L = h int ,i (29) h int ,i +1 / ,R = h int ,i +1 . (30)For piecewise linear, we set the unlimited values to h int ,i +1 / ,L, nl = h int ,i +1 / ,R, nl = (cid:26) (cid:2) ln( r i +1 /r i +1 / ) h int ,i + ln( r i +1 / /r i ) h int ,i +1 (cid:3) / ∆ ln r i +1 / (log) (cid:2) ( r i +1 − r i +1 / ) h int ,i + ( r i +1 / − r i ) h int ,i +1 (cid:3) / ∆ r i +1 / (linear)(31)and then compute the normalized slope S L = h int ,i +1 / ,L, nl h int ,i − S R = 1 − h int ,i +1 / ,R, nl h int ,i +1 . (33)We then set the final interface values to h int ,i +1 / ,L = (cid:26) h int ,i +1 / ,L, nl , | S L | ≤ lim[1 + sgn( S L )] h int ,i , | S L | > lim (34) h int ,i +1 / ,R = (cid:26) h int ,i +1 / ,R, nl , | S R | ≤ lim[1 − sgn( S R )] h int ,i +1 , | S R | > lim . (35)We adopt a fiducial value lim = 0 . VADER constructs the left andright states h int ,i +1 / ,L and h int ,i +1 / ,R from the cell-center values h int ,i usingthe piecewise parabolic method (Colella and Woodward, 1984, section 1).The reconstruction uses either r or ln r as the spatial coordinate, dependingon whether the grid is classified as linear or logarithmic. By default VADER uses piecewise linear reconstruction, as testing shows that this generally of-fers the best mix of accuracy and speed. This option is used for all the codetests presented below except where otherwise noted.10ith these approximations, the equations are fully discrete in space, andall discrete approximations are second-order accurate provided that the gridis uniform in either r or ln r . Because T depends explicitly on P , the equations to be solved are parabolicand have the form of a non-linear diffusion equation with source and sinkterms. To avoid a severe time step constraint and ensure stability, it istherefore desirable to use an implicit discretization. Let Σ ( n ) i , P ( n ) i representthe state of the system at time t n . We wish to know the state Σ ( n +1) i , P ( n +1) i at time t n +1 = t n + ∆ t . Let Θ be the time-centering parameter, such thatΘ = 0 corresponds to a forwards Euler discretization, Θ = 1 / ( n +1) i + Θ∆ t (cid:34) F ( n +1) M,i +1 / − F ( n +1) M,i − / A i − ˙Σ ( n +1)src ,i (cid:35) =Σ ( n ) i − (1 − Θ)∆ t (cid:34) F ( n ) M,i +1 / − F ( n ) M,i − / A i − ˙Σ ( n )src ,i (cid:35) (36) P ( n +1) i + Θ∆ t (cid:40) F ( n +1) P,i + − F ( n +1) P,i − A i − (cid:104) γ ( n +1) i − (cid:105) ˙ E ( n +1)int , src ,i (cid:41) = P ( n ) i − (1 − Θ)∆ t (cid:40) F ( n ) P,i + − F ( n ) P,i − A i − (cid:104) γ ( n ) i − (cid:105) ˙ E ( n )int , src ,i (cid:41) , (37)where superscript ( n ) or ( n +1) indicates whether a term is to be evaluatedusing at time t n using column density and pressure Σ ( n ) i and P ( n ) i , or at time t n +1 using column density and pressure Σ ( n +1) i and P ( n +1) i .The new time pressure fluxes F ( n +1) P,i + and F ( n +1) P,i − depend on the new spe-cific enthalpy h ( n +1) i , which in turn depends on the new internal energy E ( n +1)int ,i . For a simple equation of state with constant γ and δ = 0, we have E int = P/ ( γ − P ( r, Σ , E int ) can be inverted analyticallyto yield E int ( r, Σ , P ). However, there is no guarantee that such an analyticinversion is possible, and performing the inversion numerically could be verycomputationally costly, since it must be done in every cell for every cycle ofthe iterative method we describe below. For example, if the total pressurecontains significant contributions from both gas and radiation pressure, thenfinding the internal energy E int given a surface density and pressure wouldrequire that we solve a quartic equation. It is also possible that E int ( r, Σ , P )might not be a single-valued function of P , in which case we would face theproblem of deciding which of several possible roots is the relevant one.To avoid these problems, when γ and δ are not constant VADER evolves theinternal energy along with the column density and pressure. The evolutionof the internal energy is described by ∂E int ∂t = 1 γ − ∂P∂t + δ P Σ ∂ Σ ∂t , (38)which we discretize to E ( n +1)int ,i − E ( n )int ,i = (cid:34) Θ γ ( n +1) i − − Θ γ ( n ) i − (cid:35) (cid:104) P ( n +1) i − P ( n ) i (cid:105) + (cid:34) Θ δ ( n +1) i P ( n +1) i Σ ( n +1) i + (1 − Θ) δ ( n ) i P ( n ) i Σ ( n ) i (cid:35) (cid:104) Σ ( n +1) i − Σ ( n ) i (cid:105) . (39)This becomes our third evolution equation. By evolving the internal energyseparately we sacrifice conservation of total energy to machine precision.However, the error we make is only of order the error introduced by thediscretization of equation (38) into equation (39). We show below that ina practical example the resulting error in conservation is only marginallygreater than machine precision.For Θ = 0 the equations for the new time states are trivial, but suchan update scheme is generally unstable, or remains stable only if one obeysa time step constraint that varies as the square of the grid spacing, whichwould be prohibitively expensive in high-resolution simulations that mustrun for many viscous evolution times. For this reason, VADER supports bothΘ = 1 / / (cid:54) = 0 equations (36), (37), and (39) constitute a non-linear system,because the terms F ( n +1) E,i − / and F ( n +1) E,i +1 / that enter F ( n +1) P,i + and F ( n +1) P,i − involveproducts between P ( n +1) i +1 , P ( n +1) i , and P ( n +1) i − . If the dimensionless viscosity α , the source terms for mass ˙Σ src or internal energy ˙ E int , src , or the terms γ and δ describing the equation of state, depend on Σ or P , that represents asecond source of non-linearity. Solution when Θ (cid:54) = 0 therefore requires aniterative approach. There are numerous strategies available for solving systems of couplednon-linear equations, and our choice of method is dictated by a few con-siderations. The most common and familiar methods of solving non-linearequations are Newton’s method and its higher-dimensional generalizationsuch as Newton-Raphson iteration. However, the simplest forms of thesemethods require that we be able to compute the partial derivatives of theright-hand sides of equations (36), (37), and (39) with respect to Σ ( n ) i , P ( n ) i ,and E ( n )int ,i . Unfortunately we cannot assume that these are available, becausethe functional forms of α , ˙Σ src ˙ E int , src , γ , and δ are not known a priori , andeven if their dependence on P is known, this dependence may be tabulated,or may itself be the result of a non-trivial computation. In Section 3.3, wepresent an example application where the latter is the case. Thus we cannotassume that the derivatives of these terms are known or easily computable.While methods such as Jacobian-free Newton-Krylov (e.g. Knoll and Keyes,2004) might still be available, we instead prefer to cast the problem in termsof fixed point iteration, because it is possible to write the inner step of thisiteration in a way that is particularly simple and computationally-cheap toevaluate.Fixed point iteration is a strategy for solving non-linear equations byrecasting them as the problem of finding the fixed point of a function (Burdenand Faires, 2011, section 2.2). For example, the problem of finding the vectorof values x p that is a solution to a non-linear system of equations f ( x ) = is obviously equivalent the problem of finding a fixed point of the function g ( x ) = x − f ( x ), meaning that g ( x p ) = x p . The simplest strategy for findinga fixed point is Picard iteration, whereby one guesses an initial value x and generates a new guess by setting x = g ( x ). This procedure is then13epeated until the sequence of values x k converges to whatever tolerance isdesired. One can show that, for a well-behaved function and a starting guess x sufficiently close to the fixed point, this procedure will indeed converge.Fixed point iteration is most useful when one can make a clever choice forthe iterated function g ( x ). Note that one need not choose g ( x ) = x − f ( x ) atevery step of the iterative procedure; one only requires that, as the iterationnumber k → ∞ , g ( x ) approaches this value, or approaches any other functionthat has a fixed point when f ( x ) = 0. This means that we are free toreplace f ( x ) with a linearized function f L ( x ) that is computationally cheaperto evaluate than the true f ( x ), so long as the linearized form approaches thetrue f ( x ) as k → ∞ . This is the strategy we will adopt below.With that general discussion of fixed point iteration complete, we applythe method to our problem. Let Σ ( n ) , P ( n ) , and E ( n )int be vectors of columndensity and vertically-integrated pressure values in every cell at the start ofa time step, and let q ( n ) = ( Σ ( n ) , P ( n ) , E ( n )int ) be a combined vector describingthe full state of the simulation at time t n . We seek the vector of quantities q ( n +1) that is the solution to equations (36), (37), and (39). We can writethe formal solution as q ( n +1) = F (cid:2) q ( n ) (cid:3) , (40)where F is an operator that takes the old state q ( n ) as an argument andreturns the new state q ( n +1) .Now consider a series of guesses q ( k, ∗ ) for the true solution q ( n +1) . Wewish to generate a sequence of iterates q ( k, ∗ ) such that q ( k, ∗ ) → q ( n +1) as k → ∞ . To construct this sequence of iterates, consider a linearized versionof F , which we denote F L . The non-linear pressure equation (37) can bewritten in matrix form as M P ( n +1) = b (41)where the right hand side vector b has elements b i = P ( n ) i − (1 − Θ)∆ t (cid:40) F ( n ) P,i + − F ( n ) P,i − A i − (cid:104) γ ( n ) i − (cid:105) ˙ E ( n )int , src ,i (cid:41) + Θ∆ t (cid:104) γ ( n +1) i − (cid:105) ˙ E ( n +1)int , src ,i , (42)14nd M is a tridiagonal matrix with elements M i,i +1 = Θ ∆ tA i (cid:104) γ ( n +1) i − (cid:105) α ( n +1) i +1 (cid:40) g i +1 / r i +1 (1 − β i +1 ) (cid:34) ψ eff , i + δ ( n +1) i P ( n +1) i Σ ( n +1) i − h ( n +1) i +1 / (cid:35) + πr i +1 / v φ,i +1 / (1 − β i +1 / ) (cid:41) (43) M i,i = 1 + Θ ∆ tA i (cid:104) γ ( n +1) i − (cid:105) α ( n +1) i (cid:40) (1 − β i ) r i (cid:34) g i +1 / (cid:32) h ( n +1) i +1 / − ψ eff ,i − δ ( n +1) i P ( n +1) i Σ ( n +1) i (cid:33) + g i − / (cid:32) h ( n +1) i − / − ψ eff ,i − δ ( n +1) i P ( n +1) i Σ ( n +1) i (cid:33)(cid:35) + π (cid:2) r i +1 / v φ,i +1 / (1 − β i +1 / ) − r i − / v φ,i − / (1 − β i − / ) (cid:3)(cid:41) (44) M i,i − = Θ ∆ tA i (cid:104) γ ( n +1) i − (cid:105) α ( n +1) i − (cid:40) g i − / r i − (1 − β i − ) (cid:34) ψ eff , i + δ ( n +1) i P ( n +1) i Σ ( n +1) i − h ( n +1) i − / (cid:35) − πr i − / v φ,i − / (1 − β i − / ) (cid:41) (45) M i,j = 0 , ∀ | i − j | > . (46)The equation is non-linear because M and b both depend on the new quan-tities q ( n +1) , but we can construct a linearized version of the equation byinstead solving M L P ( † ) = b L , (47)where M L and b L differ from M and b in that all quantities that depend on q ( n +1) are replaced by identical terms evaluated with q ( k, ∗ ) . In other words,15e evaluate all terms in the matrix and in the right hand side vector usingthe previous guess at the column density and vertically-integrated pressure.Equation (47) is linear in P ( † ) . Moreover, since M L is tridiagonal, it isa particularly simple linear system, and can be solved with a number ofoperations that is linearly proportional to the number of computational cells. To complete the specification of the linearized operator F L , we linearizeequations (36) and (39) in an analogous manner, toΣ ( † ) i = Σ ( n ) i − (1 − Θ)∆ t (cid:34) F ( n ) M,i +1 / − F ( n ) M,i − / A i − ˙Σ ( n )src ,i (cid:35) − Θ∆ t (cid:34) F ( †∗ ) M,i +1 / − F ( †∗ ) M,i − / A i − ˙Σ ( †∗ )src ,i (cid:35) , (48) E ( † )int ,i = E ( n )int ,i + (cid:34) Θ γ ( k, ∗ ) i − − Θ γ ( n ) i − (cid:35) (cid:104) P ( † ) i − P ( n ) i (cid:105) + (cid:34) Θ δ ( k, ∗ ) i P ( † ) i Σ ( † ) i + (1 − Θ) δ ( n ) i P ( n ) i Σ ( n ) i (cid:35) (cid:104) Σ ( † ) i − Σ ( n ) i (cid:105) . (49)where F ( †∗ ) M,i +1 / = − g i +1 / (cid:104) α ( k, ∗ ) i +1 (1 − β i +1 ) r i +1 P ( † ) i +1 − α ( k, ∗ ) i (1 − β i ) r i P ( † ) i (cid:105) , (50)and similarly for F ( †∗ ) M,i − / . The mass source function ˙Σ ( †∗ )src is also evaluatedusing the last set of iterates for the column density, Σ ( k, ∗ ) , and the newpressure, P ( † ) , that results from solving equation (47). All the quantities onthe right hand sides of equations (48) and (49) are known, and so they canbe evaluated explicitly.Thus given a starting state q ( n ) and a guess q ( k, ∗ ) at the true solution, wecan generate a new state q ( † ) = F L (cid:2) q ( k, ∗ ) , q ( n ) (cid:3) (51)by first solving equation (47) and then solving equations (48) and (49). Sincethe linearized operator F L reduces to F if the first argument q ( k, ∗ ) = q ( n +1) , VADER ’s implementation uses the GNU Scientific Library implementation, which isbased on Cholesky decomposition (Press et al., 1992, chapter 2). See .
16t is clear that if q ( k, ∗ ) is a fixed point of F L , i.e., if q ( † ) = q ( k, ∗ ) , then q ( k, ∗ ) is also a solution to the original non-linear equation (40). Thus we haverecast the problem of solving equations (36), (37), and (39) to the problemof finding a fixed point for the linear operator F L . Recall that, if the equationof state is simple and γ and δ are constant, then we omit equation (49) anddo not evolve E int ,i . Our solution strategy can be improved further by using a convergence ac-celerator. In the Picard iteration procedure described in the previous section,one searches for the fixed point of a function g ( x ) by setting the next iterateequal to the value of the function evaluated on the previous one, i.e., by set-ting x k +1 = g ( x k ). However, this process achieves only linear convergence,meaning that, even when the guess x k is close to the true solution x p , the rateat which the residuals r k = | x k − x p | diminish obeys lim k →∞ r k +1 /r k = µ ,with µ a real number in the range (0 ,
1) (Burden and Faires, 2011, chapter10). That is, each iteration reduces the residual by a constant factor. Thegoal of a convergence accelerator is to increase the speed with which theresidual diminishes. If µ = 0, but lim k →∞ r k +1 /r qk = µ for some q > µ , then the convergence is said to be super-linear, and the goal of aconvergence accelerator is to achieve super-linearity. A side benefit of mostconvergence accelerators is to increase the radius of convergence, meaningthat the iteration procedure will still converge to the solution x p even forstarting guesses x far enough from the true solution that convergence wouldnot be achieved with Picard iteration.We choose Anderson acceleration (Anderson, 1965; Fang and Saad, 2009;Walker and Ni, 2011) as our acceleration method. The central idea of thismethod is that, rather than setting our next iterate x k +1 = g ( x k ), we in-stead set it to a weighted average of the most recent iterate and M previousones, i.e., x k +1 = (cid:80) Mj =0 ξ j g ( x k − j ). The weighting coefficients are chosen suchthat, if the function g ( x ) were linear (i.e., it had a constant Jacobian), thenthe distance between x k +1 and the exact solution would be minimized. Thiscondition can be expressed as a linear equation for the coefficients ξ j , whichcan be solved using standard matrix methods. In choosing coefficients ξ j to give the exact solution if the underlying problem were linear, Andersonacceleration functions somewhat like Newton’s method, which can jump di-17ectly to the solution in a single iteration when applied to a linear problem. While rigorous theoretical results regarding the convergence rate of Ander-son acceleration at finite M have not been derived, practical tests shows thatthe method achieves convergence rates competitive with other Newton-likemethods, which are generally super-linear (Calef et al., 2013; Willert et al.,2014).In the context of our problem, the Picard iteration method amounts tosetting q (0 , ∗ ) = q ( n ) , and then setting all subsequent iterates via q ( k +1 , ∗ ) = q ( † ) = F L (cid:2) q ( k, ∗ ) , q ( n ) (cid:3) . (52)As noted above, this strategy converges only linearly, and has a fairly smallradius of convergence, which translates into a fairly restrictive value on thetime step ∆ t . To use Anderson acceleration, rather than setting the nextiterate q ( k +1 , ∗ ) equal to F L (cid:2) q ( k, ∗ ) , q ( n ) (cid:3) , we set it to q ( k +1 , ∗ ) = M k (cid:88) j =0 ξ j F L (cid:2) q ( k − j, ∗ ) , q ( n ) (cid:3) , (53)where M k = min( k, M ), and the weight coefficients ξ j are determined byminimizing the quantity χ = (cid:88) i M k (cid:88) j =0 ( ξ j R ij ) (54)subject to the constraint (cid:80) j ξ j = 1, where R ij ≡ F L (cid:104) q ( k − j, ∗ ) i , q ( n ) i (cid:105) − q ( k − j, ∗ ) i F L (cid:104) q ( k − j, ∗ ) i , q ( n ) i (cid:105) (55)is the vector of normalized residuals in every cell after iteration k − j . For-mally, this problem is equivalent to a constrained linear least squares mini-mization of the overdetermined system R ξ = , (56) Formally, one can show that Anderson acceleration with M → ∞ is essentially equiv-alent to the generalized minimum residual (GMRES, Saad and Schultz 1986) method forsolving non-linear equations (Walker and Ni, 2011). ξ is the vector of ξ j values, is the 0 vector, and the constraintequation is · ξ = 0, where is the identity vector. This problem can besolved by a number of standard techniques. Our implementation in VADER uses QR decomposition (Press et al., 1992, chapter 2). Regardless of whether we use Picard iteration or Anderson acceleration togenerate the sequence of iterates, we repeat the calculation until the residualsatisfies max | R i, | < tol , (57)where the maximum is over all cells and all quantities in those cells after themost recent iteration, and tol is a pre-specified tolerance; VADER defaults tousing tol = 10 − , and we adopt this value for all the tests presented belowunless noted otherwise. For the system formed by the equations of mass conservation (equation7) and pressure (equation 14), we require as boundary conditions the valuesof the mass flux F M and the pressure flux F P (which depends on F M and thetorque and enthalpy fluxes F T and F E ) at the grid edge. Intuitively, we canthink of these two conditions as specifying the rate at which mass enters orexits the computational domain, and the rate at which energy enters or exitsthe computational domain.First consider the boundary condition on the mass flux. The rate at whichmass enters or exits the domain, F M , can be specified in a few ways. First,one may specify this quantity directly. However, it is often more convenientto describe the mass flux in terms of the torque. Since these two are linkedvia angular momentum conservation (equation 4), specifying the torque, orits gradient at the domain boundary, is sufficient to specify the mass flux, At present our implementation is not optimal in that we perform a QR decompositionof the residual matrix in every iteration. A faster approach would be to perform the fulldecomposition only during the first iteration, and then use QR factor updating techniquesto recompute the QR factors directly as rows are successively added to and deleted fromthe residual matrix (Daniel et al., 1976; Reichel and Gragg, 1990). However, at present noopen source implementation of the necessary techniques is available, and constructing oneis a non-trivial code development task. Since the cost of the QR decomposition is onlysignificant in problems where the viscosity and source terms are trivial to compute, andthese models are generally very quick to run in any event, we have not implemented QRfactor updating at this time. For more discussion of code performance, see Section 4. F T across the grid edge (i.e., the rate at which the applied torque doeswork on the first computational zone), or the torque in a ghost zone adjacentto the grid edge.Formally, let i = − i = N denote the indices of the ghost cells atthe inner and outer boundaries of our computational grid. The quantitiesthat we require in our update algorithm are α − P − and α N P N , since itis the combination αP that appears in the definitions of the mass, torque,and enthalpy fluxes. Without loss of generality we can set α − = α and α N = α N − , since only the combination αP matters. With this choice, weset P − and P N as follows:1. If the mass flux F M, − / or F M,N − / across the grid edge is specified,then from equation (20) we have P − = (1 − β ) r (1 − β − ) r − P + F M, − / g − / α − (1 − β − ) r − (58) P N = (1 − β N − ) r N − (1 − β N ) r N P N − − F M,N − / g N − / α N (1 − β N ) r N . (59)where r − = e − ∆ ln r r for a logarithmic grid, or r − = r − ∆ r for alinear one, and similarly for r N .2. If the torque flux F T, − / or F T,N − / is specified, then from equation(22) we have P − = − P + F T, − / πr − / v φ, − / (1 − β − / ) α − (60) P N = − P N − + F T,N − / πr N − / v φ,N − / (1 − β N − / ) α N . (61)3. If the torque T − or T N in the first ghost zone is specified, then fromequation (5) we have P − = − T − πr − (1 − β − ) α − (62) P N = − T N πr N (1 − β N ) α N . (63)20t should be noted that the inner and outer boundary conditions need notbe specified in the same way, e.g., one can specify the mass flux indirectlyby giving the torque at the inner boundary, and set the mass flux directly atthe outer boundary.Now consider the second boundary condition we require. Intuitively thisis the rate of energy transport into and out of the domain, but since we areworking in terms of a pressure equation, it more convenient to think in termsof the flux of enthalpy F E into or out of the domain. From equation (23), thisdepends on both the mass flux, and hence αP , and on the specific internalenthalpy h int . Thus we specify the final part of the boundary conditions bysetting h int , − and h int ,N , the values of internal enthalpy in the ghost cells.Intuitively, if we know the mass flux into the domain from the first boundarycondition, and we have now specified the enthalpy of the material that isbeing advected in, we can compute the rate at which advection is bringingenergy into the computational domain. The enthalpies we require can ei-ther be specified directly, or specified in terms of the gradient of enthalpyacross the grid boundary. Note that if mass is always flowing off the gridat either the inner or outer boundary, and if one adopts piecewise constantinterpolation, then the boundary value chosen for h int will not affect the solu-tion. Also note that, because we describe the boundary condition by setting h int , − rather than by setting F E, − / directly, the boundary specification isindependent of the method chosen to reconstruct the enthalpies at cell edges.A final subtlety in choosing boundary conditions is that these conditionsmust be applied to both the old and new times, i.e., they must apply toboth P ( n ) and P ( n +1) . The iterative solver will ensure that this holds tothe accuracy of the iterative solve regardless, but by assigning appropriatevalues to the parts of the right hand side vector b and matrix M that affectthe ghost cells (i.e., elements M − , − , M − , , M N,N − , M N,N , b − and b N ) onecan enforce the boundary conditions to machine precision at both the oldand new times, and in the process speed convergence of the iterative solutionstep. The boundary conditions take the form of a relationship between theghost cells and the adjacent real cells, as specified by equations (58) – (63), The enthalpy flux F E also depends on the effective gravitational potential ψ eff ,i acrossthe grid boundary, but we will assume that this is known from the specification of therotation curve. In general it is also possible to have boundary conditions such that values in theghost cells depend on the values in real cells that are not immediately adjacent to the P − = q − P + p − (64) P N = q N P N − + p N , (65)with the constants q and p depending on how the boundary condition isspecified. For a given choice of q and p , the matrix and vector elementsrequired to enforce the boundary conditions are M − , − = M N,N = 1 (66) M − , = − q − (67) M N,N − = − q N (68) b − = p − (69) b N = p N . (70) The final piece required to complete the update algorithm is a recipe tochoose the time step ∆ t . Using either Θ = 1 / / VADER controls the time stepbased on the rate of change computed in the previous time step. If q ( n − was the state of the system at time t n − , q ( n ) is the state at time t n , and∆ t n − was the previous time step, then the next time step ∆ t n = t n +1 − t n is determined by∆ t n = C min i =0 , ,...N − (cid:32)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) q ( n − i q ( n ) i − q ( n − i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:33) ∆ t n − . (71)where the minimum is over all cells, and C is a dimensionless constant forwhich we choose a fiducial value of 0.1. To initialize the calculation, VADER ghost zones. In this case the best approach is to treat the boundary condition as simplya specified value of P − or P N , without attempting to enforce this relationship at boththe old and new times by setting elements of M . Including these constraints in M wouldrender M no longer tridiagonal. Since this would increase the cost of solving the linearsystem considerably, it is more efficient to include only the relationship between the ghostcells and their nearest neighbors in M , and enforce any other relationships between theghost and real cells by iterating the system to convergence. − r − / /v φ, − / , during which the code com-putes a new state and uses it to determine the time step, but does not actuallyupdate the state.While the update algorithm is stable regardless of the choice of time step,the same is not necessarily true of the iterative solution procedure describedin Section 2.3, which may not converge if ∆ t is too large. The maximum timestep for which convergence occurs will in general depend on the boundaryconditions and the functional forms of α , γ , δ , ˙Σ src , and ˙ E src . To ensure thestability of the overall calculation, VADER allows a user-specified maximumnumber of iterations used in the implicit solve. If convergence is not reachedwithin this number of iterations, or if the iteration diverges entirely,
VADER stops iterating and tries to advance again with a factor of 2 smaller timestep. This ratcheting is applied recursively as necessary until convergence isachieved.
We have implemented the algorithm described above in the Viscous Ac-cretion Disk Evolution Resource (
VADER ), and made the code publicly avail-able from https://bitbucket.org/krumholz/vader/ under the terms ofthe GNU Public License. The core code is written in c, with templates pro-vided for user-supplied implementations of functions characterizing the innerand outer boundary conditions, the dimensionless viscosity, the equation ofstate terms, and the rates of mass and energy gain or loss per unit area. Eachof these quantities can also be specified simply by a numerical value, so thatusers who do not wish to use some particular aspect of the code (e.g. thegeneral equation of state) need not implement functions describing it.
VADER also provides routines to construct rotation curves suitable for computationaluse from tabulated v φ ( r ) values; see Appendix A for details.In addition to the core C routines, VADER includes a set of Python wrapperroutines that allow simulations to be run from a Python program. ThePython wrappers provide high-level routines for processing simulation inputand outputs, and controlling execution.The
VADER repository includes an extensive User’s Guide that fully doc-uments all functions. 23 . Accuracy and Convergence Tests
In this section we present tests of the accuracy and convergence charac-teristics of the code, comparing against a variety of analytic solutions, anddemonstrating a number of the code’s capabilities, including linear and log-arithmic grids, a range of rotation curves, and arbitrary functional forms forenergy gain and loss, boundary conditions, and viscosity. The code requiredto run all the tests described in this section is included in the
VADER bitbucketrepository.One subtlety that occurs in all these comparisons is treatment of theboundary conditions. The true nature of the boundary layers of accretiondisks is subject to significant theoretical uncertainty (e.g., Papaloizou andStanley, 1986; Popham and Narayan, 1992; Popham et al., 1993), but mostknown analytic solutions for viscous disks generally prescribe boundary con-ditions by specifying that the disk extends all the way from r = 0 to r = ∞ ,and by demanding regularity as r → r → ∞ . (For an exceptionsee Tanaka (2011), who derives analytic solutions for disks with power lawviscosities and finite inner radii.) For obvious reasons in a numerical compu-tation (and in nature) it is necessary to truncate the disk at finite values of r ,and thus some care is required to match the boundary conditions used in thenumerical computation to those assumed in the analytic solutions. For mostreasonable choices of boundary condition, the choice will affect the resultsonly in the vicinity of the boundaries, but for the purpose of making quantita-tive comparisons one must be careful to match boundary conditions betweenthe analytic and numerical solutions to ensure that differences caused bydifferent boundary conditions do not dominate the error budget. The first test is a comparison to the similarity solution derived by Lynden-Bell and Pringle (1974). The solution is for Keplerian rotation ( β = − / ν that varies with radius as ν = ν ( r/R ), giving α = (cid:18) ν v φ R (cid:19) Σ P . (72)Note that this choice of α renders the torque a function of Σ only, so thatthe transports of mass and energy are decoupled and the rates of transportdo not depend on P . With this rotation curve and viscosity, the evolution24quations admit a similarity solutionΣ = Σ e − x/T xT / , (73)where Σ = ˙ M / (3 πν ), x = r/R , T = t/t s , t s = R / ν , and ˙ M is themass accretion rate reaching the origin at time T = 1. To check VADER ’sability to reproduce this solution, we simulate a computational domain from r/R = 0 . −
20, using a uniform logarithmic grid. We initialize the grid usingthe analytic solution at time T = 1, and specify the boundary conditions bysetting the torque in the ghost zones equal to the values for the similaritysolution, T = − ˙ M v φ R xT / e − x/T . (74)The equation of state is a simple one, with γ = 1 + 10 − and δ = 0, and theinitial value of P/ Σ is set to a constant. The value of γ is chosen so that P/ Σ remains nearly constant.Figure 1 shows a comparison between the analytic solution and the resultsof a
VADER simulation performed with a resolution of 512 cells, piecewise-linear reconstruction, a tolerance C = 0 . − numerical − Σ analytic Σ analytic , (75)is generally of order 10 − − − , with a maximum of about 5 × − at time t/t s = 2. Note that a small amount of grid-scale oscillation is visible in theerror at small radii, although the overall error is still ∼ − . Oscillationsof this sort are a common feature of calculations using Crank-Nicolson timecentering, and can be eliminated by use of the backwards Euler method, atthe price of lower overall accuracy. The magenta line in Figure 1 shows anotherwise identical calculation using backwards Euler. The grid noise in theCrank-Nicolson result can also be greatly reduced by using a more stringenttolerance in the iterative solve. Using tol = 10 − instead of the default valueof 10 − , renders the oscillation invisibly small. This is shown by the cyanline in Figure 1.To determine how the accuracy of the solution depends upon spatial res-olution, we repeat the calculation using Crank-Nicolson time centering at25 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 Σ / Σ Analytic t/t s =1 t/t s =2 t/t s =3 t/t s =4 Numerical t/t s =1 t/t s =2 t/t s =3 t/t s =4 -1 r/R -7 -6 -5 -4 -3 -2 -1 | E rr o r | t/t s =2 t/t s =3 t/t s =4 t/t s =4 , BE method t/t s =4 , tol = − Figure 1: Comparison between the analytic solution for the self-similar evolution of aviscous disk (equation 73) and a
VADER simulation. The upper panel shows the normalizedgas surface density Σ / Σ versus normalized radius r/R at several times for the analyticsolution (solid lines) and the simulation (circles; only every eigth data point plotted toavoid confusion). The values plotted for t/t s = 1 are the initial conditions in the simu-lation. The lower panel shows the absolute value of the error, defined per equation (75).The red, green, and blue lines show the same calculation as in the upper panel. Themagenta line shows the result at t/t s = 4 using a backwards Euler method instead of aCrank-Nicolson one, while the cyan line shows a result using a Crank-Nicolson methodbut with a tighter error tolerance of 10 − instead of 10 − in the iterative solution step.
26 range of resolutions from N = 64 to 2048 in factor of 2 steps. We usetol = 10 − for this test to ensure that the error is determined by the spatialresolution, and not the choice of error tolerance in the iterative solver. Figure2 shows the absolute value of the error versus position at a time t/t s = 2 forthese calculations, and Figure 3 shows the L error in the solution at thistime versus resolution, where the L (Lebesgue) norm of the error has itsusual definition L error = 1 π Σ R (cid:88) i A i | Σ numerical ,i − Σ exact ,i | . (76)As expected, the error declines with resolution, with the L error declining as N − . A least-squares fit of the logarithm of the L error versus the logarithmof resolution gives a slope of − .
0. This confirms that the code is second-order accurate in space, as expected. We do not perform a similar test ofthe time accuracy, because the actual simulation time step is not controlledby a single parameter. It depends on the time step tolerance C , but alsoon the error tolerance tol used in the iterative solver, and on the maximumnumber of iterations the solver is allowed to perform before giving up andtrying again with a reduced time step. The second test we present is a comparison of
VADER to the analytic so-lution for the evolution of an initially-singular ring of material with constantkinematic viscosity (Pringle, 1981). Consider a disk whose initial densitydistribution is concentrated at a single radius r = R , such thatΣ( r, t = 0) = m δ ( r − R )2 πR . (77)The material orbits in a Keplerian potential, and has a constant kinematicviscosity ν , corresponding to α = ν (cid:18) Σ P (cid:19) (cid:16) v φ r (cid:17) . (78)As in the previous test problem, this choice of α renders the torque a functionof Σ only. With this choice of α , and subject to the boundary conditions that27 -1 r/R -8 -7 -6 -5 -4 -3 -2 -1 | E rr o r | N =64 N =128 N =256 N =512 N =1024 N =2048 Figure 2: Same as the lower panel of Figure 1, but now showing the absolute value of theerror versus position at time t/t s = 2 for a series of computations with different numbersof cells N . All the calculations shown use Crank-Nicolson time centering and iterativesolver tolerance tol = 10 − . N -7 -6 -5 -4 -3 -2 L E rr o r Numerical N − Figure 3: L error (equation 76) versus resolution N at time t/t s = 2 for the self-similardisk test results shown in Figure 2 (blue lines and points). The black dashed line shows aslope of −
29 remain finite and
T → r → the system has the exact solutionΣ exact = Σ x / τ e − (1+ x ) /τ I / (2 x/τ ) , (79)where x = r/R , τ = t/t s , Σ = M /πR , and I / is the modified Besselfunction of the first kind of order 1 /
4. Here t s = r / ν is the characteristicviscous evolution time.Due to the singular initial condition, this is a far more challenging testof the algorithm than the similarity solution discussed in the previous sec-tion. To test VADER ’s performance on this problem, we simulate a disk ona uniform linear grid extending from 0 . R − R , initialized such that thecell containing the radius R has a column density Σ = Σ init ≡ M /A , where A is the area of the cell. All other cells have column densities Σ = Σ init /χ with χ = 10 . The initial vertically-integrated pressure is set so that P/ Σis constant, and the simulation uses an equation of state with γ = 1 + 10 − , δ = 0, but neither of these choices affects the evolution of the gas surfacedensity. We set the torques at the inner and outer boundaries equal to theiranalytic values (subject to the density floor), T = − πrνv φ Σ , (80)where Σ = max(Σ exact , Σ init /χ ). Note that the boundary torque is thereforetime-dependent. Figure 4 shows the results of the test for a simulation with4096 cells, which required 20 −
30 seconds of run time to complete. As theplot shows,
VADER reproduces the analytic result very accurately. The error,defined to account for the effects of the density floor asError = Σ numerical − Σ exact − Σ init /χ Σ exact + Σ init /χ , (81)is of order 10 − at very early times when the ring is poorly resolved, anddrops to ∼ − or less at late times. The third test we present is against the analytic solution for gravita-tional instability-dominated disks first derived in a special case by Bertin Because the equations for Σ and P are decoupled, only two boundary conditions arerequired to specify the solution for Σ. -6 -5 -4 -3 -2 -1 Σ / Σ Analytic t/t =0 . t/t =0 . t/t =0 . t/t =0 . Numerical t/t =0 . t/t =0 . t/t =0 . t/t =0 . r/R -6 -5 -4 -3 -2 E rr o r Figure 4: Results of a simulation of an initially-singular ring test. The upper panel showsthe exact analytic solution (solid lines, equation 79) and the numerical results producedby
VADER (circles; only every 64th point shown, for clarity) as a function of position r/R at times t/t = 0 . VADER ’s ability to handlea much more complex problem than the previous two tests. In this problemthe rotation curve is not Keplerian, radiative cooling is an integral part ofthe problem, and the viscosity is not given by a simple analytic formula butinstead is derived from an underlying set of equations to be solved at eachtime step.Consider a pure gas disk (i.e., one with no stellar component) of surfacedensity Σ, within which support against self-gravity is dominated by a highlysupersonic velocity dispersion σ , such that the vertically-integrated pressure P = Σ σ , γ = 5 /
3, and δ = 0. The stability of the disk against axisymmetricgravitational instabilities is described by the Toomre (1964) Q parameter, Q = κσ/πG Σ, where κ = [2( β + 1)] / v φ /r is the epicyclic frequency. Theturbulence in the disk decays following˙ E int , src = − η Σ σ v φ r , (82)where η = 3 / Q is relatedto the torque implicitly via τ (cid:48)(cid:48) + h τ (cid:48) + h τ = H, (83)where τ = T / [ ˙ M R v φ ( R ) R ] is the dimensionless torque, R is a chosen scaleradius, v φ ( R ) and ˙ M R are the rotation curve speed and initial (inward) ac-cretion rate at that radius, and the coefficients appearing in the equation Formally Toomre’s Q applies only to a disk where σ is the thermal velocity dispersion,but a generalized Q applies to gas where the non-thermal velocity dispersion is muchgreater than the sound speed (e.g., Elmegreen, 2011). h = ( β − u x s (84) h = − β + 1) xs (cid:48) + 2 s ( β + β + xβ (cid:48) )2( β + 1) sx (85) H = (cid:115) ( β + 1) π χ (cid:18) πηu − d ln Qd ln T (cid:19) su Qx (86)where x = r/R , u = v φ /v φ ( R ), s = σ/v φ ( R ), χ = G ˙ M R /v φ ( R ) , T = t/t orb , t orb = 2 πR/v φ ( R ), and primes indicate differentiation with respect to x .If d ln Q/d ln T = 0 for Q = 1, and β is constant, then the combinedsystem formed by equation (83) and the equations of mass and energy con-servation (1) and (2) admit a steady-state solution. For a flat rotation curve,i.e., one with v φ constant and β = 0, this is T = − r ˙ M R v φ ( R ) (87)Σ = v φ πGr (cid:32) G ˙ M R η (cid:33) / (88) σ = 1 √ (cid:32) G ˙ M R η (cid:33) / . (89)If d ln Q/d ln T ≤ Q >
Q >
VADER ’s ability to solve this problem, we perform a series oftests. In all of these simulations we obtain the viscous torque T and thus thedimensionless viscosity α required by VADER by solving a discretized versionof equation (83) on the grid with d ln Qd ln T = ux min( e − /Q − e − , , (90)so that the disk returns to Q = 1 on a timescale comparable to the localorbital time. This particular functional form is not particularly physicallymotivated, and is chosen simply to ensure that d ln Q/d ln T goes smoothly33o 0 as Q → τ = − x at the inner boundary and τ = − β − Q >
1, wesuppress gravitational instability-driven transport by reducing the torque bya factor e − Q − . The VADER simulation uses an outer boundary conditionwith a fixed mass flux ˙ M , and an inner boundary condition whereby the fixedtorque is given by τ = − xe − Q − , where Q is evaluated in the first gridzone. All simulations use a uniform logarithmic grid of 512 cells that goesfrom 0 . R − R , and piecewise constant enthalpy advection.Figure 5 shows the results of a simulation where the system is initializedto the analytic steady-state solution and allowed to evolve for 4 outer orbitaltimes, corresponding to 400 orbital times at the inner edge of the disk. Asthe figure shows, VADER successfully maintains the steady state. The totaltime required to run this computation was 20 −
30 seconds.Although no analytic solutions are known for systems that start awayfrom equilibrium, a more stringent test is to start the code away from theanalytic solution and verify that it approaches the steady state in a physicallyreasonable manner. Figure 6 shows the evolution for a simulation that beginsthe same surface density as the steady state solution (equation 88), but avelocity dispersion a factor of 2 smaller, and thus at Q = 0 .
5. The enthalpyat the outer boundary condition is also a factor of 2 below the steady-statesolution value. As the Figure shows, the system rapidly evolves from theinside out to Q = 1 due to an increase in the velocity dispersion drivenby viscous transport. After ∼ ∼ Q = 1,but with the surface density and velocity dispersion both increased by a fac-tor of 2 relative to the steady state solution. The outer enthalpy boundarycondition is equal to that of the analytic solution. Again, the system evolvestoward the steady-state solution. The time required for this evolution is ∼ ∼ − Σ / Σ ( R ) Steady stateSimulation, T =0 . Simulation, T =4 . -2 -1 σ / v φ -2 -1 r/R Q Figure 5: Results of a simulation of the gravitational instability-dominated disk test.The three panels shows the gas surface density Σ normalized to the steady-state value at R , the middle panel shows the velocity dispersion σ normalized to v φ , and the bottompanel shows Toomre Q . The black dashed line shows the analytic steady state solution(equations 88 and 89), the blue line shows the simulation initial condition, and the greenline shows the simulation after T = 4 outer orbits. Note that the blue line is completelyhidden by the green line, as it should be since we are testing the ability of the code tomaintain the correct analytic steady state. Σ / Σ ( R ) Steady stateSimulation, T =0 . Simulation, T =0 . Simulation, T =0 . Simulation, T =1 . -2 -1 σ / v φ -2 -1 r/R Q Figure 6: Same as Figure 5, but now showing a simulation that starts out of equilibriumwith Q = 0 . Σ / Σ ( R ) Steady stateSimulation, T =0 . Simulation, T =0 . Simulation, T =10 . Simulation, T =20 . -2 -1 σ / v φ -2 -1 r/R Q Figure 7: Same as Figure 5, but now showing a simulation that starts out of equilibriumwith Q = 1 but a surface density and velocity dispersion that are both double the steady-state value. .4. Singular Ring with a Complex Equation of State Our fourth and final test focuses on
VADER ’s ability to handle complexequations of state. As is the case for a non-equilibrium gravitational instability-dominated disk, no analytic solutions are known for this case, but we cannevertheless verify that the code gives a physically realistic solution, andthat it shows good conservation properties. We therefore choose to repeatthe singular ring test, but using parameters such that the ring encountersboth gas pressure-dominated and radiation pressure-dominated regimes.For simplicity consider a disk where the shapes of the vertical density andtemperature distributions are fixed, so that the density and temperature canbe separated in r and z , i.e., ρ ( r, z ) = ρ ( r ) a ( ζ ) (91) T ( r, z ) = T ( r ) t ( ζ ) (92)where ζ = z/z is a dimensionless height, z is an arbitrary vertical scalefactor, and the vertical distribution functions a ( ζ ) and t ( ζ ) are both normal-ized such that (cid:82) a ( ζ ) dζ = (cid:82) t ( ζ ) dζ = 1. The vertically-integrated columndensity and gas and radiation pressures areΣ = ρz (93) P g = k B µm H Σ T (cid:90) a ( ζ ) t ( ζ ) dζ ≡ k B µm H Σ T eff (94) P r = 13 aT z (cid:90) t ( ζ ) dζ ≡ f aT z , (95)where µ is the mean molecular weight, and we have defined T eff = T (cid:90) a ( ζ ) t ( ζ ) dζ (96) f = (cid:82) t ( ζ ) dζ (cid:2)(cid:82) a ( ζ ) t ( ζ ) dζ (cid:3) . (97)The corresponding vertically-integrated internal energies are E int ,g = P g γ g − E int ,r = 3 P r , (99)38here γ g is the gas adiabatic index. The total pressure and internal energyare simply the sums of the gas and radiation components, i.e., P = P g + P r and E int = E int ,g + E int ,r . From these definitions, with a bit of algebra onecan show that the equation of state parameters γ and δ can be written interms of the total pressure and internal energy as γ = (16 − γ g ) P + (16 − γ g ) E int P + (13 − γ g ) E int (100) δ = 4(3 P − E int ) [( γ g − E int − P ] P [3( γ g − E int + (3 γ g − P ] . (101)Since the problem is not dimensionless once a real equation of state isadded, we adopt the following dimensional parameters. The simulation do-main has inner radius r − / = 1 . × cm, outer radius r N − / = 1 . × cm, and a rotation curve corresponding to Keplerian motion about a cen-tral object of mass M = 3 M (cid:12) . The initial ring of material is located at R = 7 . × cm, its mass and effective temperature are M = 1 . × − M (cid:12) and T eff = 10 K, its adiabatic index γ g = 5 /
3, and the scale heightparameter f z = 7 . × cm. The kinematic viscosity ν = 1 . × cm s − , which gives a characteristic evolution time t s = 10 yr. With thesechoices, the ring has P g (cid:29) P r in its interior, but as it expands, its edgesheat up to the point where P r (cid:29) P g . The simulations use a linear grid of4096 cells, backwards Euler updating, and piecewise constant interpolationof enthalpy; the latter two choices are made to suppress numerical ringing atthe discontinuous expansion front. We also use a time step tolerance C = 1 . . ∼ ∼ γ = 5 / δ = 0, and a case including radiation pressure. (Thecolumn density distributions are nearly identical to those shown in Figure4, as is expected since the choice of viscosity for this problem is such thatthe column density evolution does not depend on the pressure.) As shownin the Figures, in both cases the spreading ring has sharp rises in pressureand temperature at its low-density edges, where the ring encounters near-vacuum. However, in the run including radiation pressure the pressure andtemperature in these spikes are greatly reduced. In the run without radia-tion pressure, the temperature rises to unphysical values because the pressure39ust be carried entirely by the gas, and the gas pressure is linear in tem-perature and inversely proportional to surface density. In the run includingradiation, the pressure at the ring leading and trailing edges is carried byradiation instead, and the temperature rise is vastly reduced because the ra-diation pressure rises as the fourth power of temperature, and does not scalewith gas surface density. There is also a small spike in the pressure at theoriginal ring location in the case with radiation pressure. It is not clear ifthis is physically real, or if the spike is a purely numerical effect resultingfrom the unresolved, singular initial condition.This test also enables us to check the level of conservation of energy in thecomputation with the complex equation of state. We simulate a time intervalfrom t/t s = 0 to 0 . γ we find that the maximum change in totalenergy in the computational domain, after accounting for energy transmittedacross the grid boundary by advection or torques, is 9 . × − of the initialenergy, and that the mean difference between the initial energy and theenergy at later times is 4 . × − of the initial energy. This is consistent withour expectations that the algorithm should conserve total energy to machineprecision. In contrast, the run including radiation pressure has maximumand mean energy conservation errors of 3 . × − and 1 . × − of theinitial energy. Thus, while conservation is not quite at machine precision,the loss of precision is only ∼
4. Performance
The tests presented in the previous section demonstrate that
VADER andthe algorithms on which it is based provide correct solutions to a numberof problems, and give a rough indication of the performance of the code.Here we investigate the performance of the code in much more detail. Weare particularly interested in the performance of the implicit solver and theAnderson acceleration code, because, while Anderson acceleration acceleratesconvergence and reduces the number of iterations requires in the implicitsolver, it also requires a linear least squares solve that increases the costper iteration. The tradeoff between these two is almost certainly problem-dependent, and may also be processor- and compiler-dependent, but thetests we describe here can serve as a guide for users in selecting appropriateparameters for their own problems. All the tests we discuss in this section40 P / f z [ d y n c m − ] Gas only ( γ =5 / ) t/t =0 . t/t =0 . t/t =0 . r/R P / f z [ d y n c m − ] Gas plus radiation
TotalGasRadiation
Figure 8: Vertically-integrated pressure P divided by the scale height parameter f z atseveral times in the singular ring test. The top panel shows a simulation omitting radiationpressure ( γ = 5 / δ = 0), while the bottom panel shows an otherwise identical simulationwith a complex equation of state including radiation pressure. In the gas plus radiationrun, we show both total pressure (solid lines) and the pressures due to gas (dashed lines)and radiation (dotted lines) alone. T e ff [ K ] Gas only ( γ =5 / ) t/t =0 . t/t =0 . t/t =0 . r/R T e ff [ K ] Gas plus radiation
Figure 9: Same as Figure 8, but now showing effective temperature T eff versus positionin the singular ring test without (top) and with (bottom) radiation pressure. Notice thedifference in scales between the top and bottom panels. VADER was compiled using gcc-4.8 with optimizationlevel -O3 , while the GNU Scientific Library was built using its default options.We obtain code timing using the C clock() function.We first verify that Anderson acceleration does, as expected, lead to muchmore rapid convergence of the iterative solver. To test this, we run each of ourfour test problems described in Section 3 – the self-similar disk, the singularring, the gravitational instability-dominated disk, and the radiation pressurering – for one time step, starting from the initial conditions as described inthe previous Section. We use time steps of 10 − . t s , 10 − t s , 10 − . t orb , and10 − . t s , respectively, and we test both Crank-Nicolson and backwards Eulerupdating. We set the tolerance on the iterative solver to 10 − , and allowa maximum of 100 iterations. Figure 10 shows how the residual changesversus number of iterations for each of these runs, and Table 1 shows thenumber of iterations and the wall clock time required to converge. Theperformance of the algorithm is in line with our expectations. We also seethat the Crank-Nicolson method almost always produces faster convergencethan the backwards Euler method, which is not surprising due to its higherorder of accuracy.Note that, at lower orders of Anderson acceleration, for some problems theiteration diverges rather than converging, until eventually the code producesnon-numerical values ( Inf or NaN ), at which point the solver halts iteration.In a full simulation, such cases of divergence are treated exactly like caseswhere the solver fails to converge within the prescribed maximum number ofiterations, i.e., the time step is attempted again using a reduced value of ∆ t (see Section 2.5).While this test shows that Anderson acceleration does speed convergencein terms of number of iterations, and does allow larger time steps, it does notprove that the extra computational cost per iteration is worthwhile. Indeed,carefully examining the timing results given in Table 1, we see that thewall clock time is by no means a monotonically decreasing function of M ,even in cases where the number of iterations is. To evaluate this question,we next run each of our test simulations for 1000 time steps or until thesimulation completes, whichever comes first, and measure the total executiontime. For this test we return the iterative solver tolerance to its default valueof 10 − , the maximum number of iterations to its default of 40, and allowthe time step to be set by the normal VADER procedure. We use a time steprestriction C = 0 . -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 R e s i d u a l Self-similar Ring Crank-Nicolson M =0 M =1 M =2 M =4 M =8 M =16 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 R e s i d u a l GI Disk 0 20 40 60 80 100Iteration numberRad. Ring Backwards Euler M =0 M =1 M =2 M =4 M =8 M =16 Figure 10: Maximum normalized residual max | R ,i | (equation 57) remaining after N iterations in the four test problems described in Section 4. Solid lines show updates usingthe Crank-Nicolson method, and dashed lines using the backwards Euler method. Colorsindicate the order M of Anderson acceleration used, with M = 0 corresponding to noacceleration (standard Picard iteration). elf-similar Ring GI Disk Rad. ring M CN BE CN BE CN BE CN BE N iter . . . . . .
22 – 94 – – – 60 –1 . . . . . .
21 – 92 – – – 26 252 . . . . . .
16 51 44 – 72 100 26 244 . . . . . .
14 31 44 72 27 89 26 258 . . . . . .
13 23 39 50 22 31 28 3116 . . . . . .
13 20 31 48 21 26 27 36Time [ms] 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Table 1: Number of iterations and total wall clock time required for convergence duringone step of the self-similar disk, singular ring, gravitational instability-dominated disk,and singular ring with radiation pressure tests, using Crank-Nicolson and backwards Eulerupdating methods. The quantity M is the Anderson acceleration parameter, with M = 0indicating no acceleration (standard Picard iteration). Blank entries indicate that thesolver failed to converge within 100 iterations. and C = 1 . M and actually harm-ful at larger M . This is because the reduction in number of iterations ismore than offset by the increased cost per iteration. On the other hand, forthe gravitational instability-dominated disk and ring with radiation pressureproblems, Anderson acceleration with M of a few provides very significantgains in performance, reducing the computational cost by a factor of ∼ M . In contrast, the gravitationalinstability-dominated disk and radiation pressure ring problems have muchhigher computational costs per update. For the former, the cost is high be-cause there is a source term and because computing the viscosity requires45olving a tridiagonal matrix equation at every iteration; this is indicated bythe green region in Figure 11. For the latter there are no source terms, butbecause the problem uses a complex equation of state, there is additionalcomputational work associated with updating the internal energy equation(which is included in the red region in Figure 11).The implications of this analysis are that the choice of optimal Andersonacceleration parameter is likely to depend on the nature of the model beingused to generate the viscosity, source terms, boundary conditions, and equa-tion of state. If these are given by simple analytic formulae, then Andersonacceleration is probably of very limited use. The more computationally com-plex they are to evaluate, however, the greater the advantage that one gainsby using Anderson acceleration to reduce the number of iterations required.
5. Summary and Future Prospects
Thin, axisymmetric accretion disks where the transport of mass and an-gular momentum is approximated as being due to viscosity represent animportant class of models in theoretical astrophysics. They are widely usedin situations where full two- or three-dimensional simulations would be pro-hibitively expensive, either due to the number of orbits that would have tobe simulated, or because of the dynamic range in spatial resolution thatwould be required. We have developed a new, extremely flexible and gen-eral method for simulating the evolution of such models. Our discretizationof the equations is conservative to machine precision, and allows completefreedom in the specification of rotation curves, equations of state, formsof the viscosity, boundary conditions, and sources and sinks of energy andmass. The core of our method is an unconditionally-stable update strat-egy that uses accelerated fixed point iteration to achieve rapid convergence.We show that this technique allows relatively large time steps, and that itsignificantly reduces the number of implicit iterations required to advancethe simulation a specified time. In practical tests, even at resolutions ashigh as 4096 cells, and using complex equations of state or disk α param-eters that must themselves be computed iteratively, the code is can evolvea disk for many viscous evolution timescales in no more than a few min-utes of computation time on a single CPU. In tests with simple equationsof state and fixed α , computational times are only a few seconds. We haveimplemented our algorithm in a new open source code called the ViscousAccretion Disk Evolution Resource VADER , which can be downloaded from46 .00.51.01.52.02.53.0 C o m p u t a t i o n a l c o s t Self-similar
TotalProblem-specificAndersonAdvance
Ring0 1 2 4 8 M C o m p u t a t i o n a l c o s t GI Disk 0 1 2 4 8 16 M Rad. Ring
Figure 11: Wall clock time required per unit simulation time advanced versus Andersonacceleration parameter M in the four test problems, normalized so that M = 0 correspondsto unity. Thus values < > α ), source terms ( ˙Σ src , ˙ E int , src ), equation of state terms ( γ , δ ), and boundaryconditions are shown in green. The Anderson acceleration step is shown in blue; for M = 0 this is simply the cost of copying temporary arrays. The remainder of the timestep advance procedure is shown in red. The dashed horizontal lines indicate values of0 .
25, 0 .
5, 1 .
0, and 2 .
0. Note the change in y -axis range between the top two and bottomtwo panels. The plots shown are for Crank-Nicolson; the results for backwards Euler arequalitatively similar. ttps://bitbucket.org/krumholz/vader/ . The code is designed for mod-ularity, so that users can easily implement their own models for viscosity andsimilar parameters that control disk evolution.While the number of potential applications for VADER is large, we end thisdiscussion by highlighting a few possible examples. One, already underway(Krumholz & Kruijssen, 2015, in preparation), is to model the long-termbehavior of the gas around the central black in the Milky Way and othergalaxies. Observations suggest that the gas accumulates over long timescalesbefore undergoing periodic starburst events (e.g. Kruijssen et al., 2015), andthis process can be modeled as gas undergoing slow viscous accretion beforeaccumulating to the point where it becomes gravitationally unstable andundergoes a starburst. The customization required for this project consistsmostly of implementing a custom rotation curve, viscosity representing theinstabilities that likely drive the gas migration, and testing for the onset ofgravitational instability.In the context of planet formation, viscous evolution models have beenused to study the long-term interaction of planets with the disks out ofwhich they form, and the migration of the planets through viscous disks(e.g., Lyra et al., 2010; Horn et al., 2012).
VADER is well-suited for thisapplication, and could be customized to it simply by implementing coolingand viscosity prescriptions appropriate to a protoplanetary disk, by couplingthe state of the disk found by
VADER to a calculation of the torques on planets,and perhaps by modifying the viscosity prescription to incorporate the back-reaction of disk-planet torques on the transport of gas through the disk. Alsoin this context,
VADER could be used to simulate the photoevaporation of disksaround young stars by high-energy radiation, either from the central star orfrom an external source (e.g. Adams et al., 2004; Gorti and Hollenbach, 2009;Gorti et al., 2009). In this case one could treat the effects of photoevaporationby adding mass and energy source terms to represent the rates of mass lossand heating driven by ultraviolet and X-ray photons striking the disk.For accretion disks around compact objects, a number of authors haveused viscous disk models to study variability and flaring on timescales as-sociated with viscous evolution (e.g., Cambier and Smith, 2013).
VADER iswell-suited to this problem too, and could be customized to it by adding inprescriptions for viscosity and radiative heating and cooling, and by post-processing the
VADER models to predict observable X-ray fluxes. One couldalso add mass and energy source terms representing the exchange of massand energy with a hot corona. 48his list of potential applications is certainly not exhaustive, but itsbreadth and diversity should make clear that a general-purpose viscous diskevolution code is a tool of wide applicability. The availability of such a codeshould reduce the need for every modeler to develop his or her own approachto a standard problem.
Acknowledgements
We thank D. Kruijssen for stimulating discussions that helped inspirethis work, and the authors and maintainers of the GNU Scientific Library forproviding a valuable resource. We thank W. Lyra and T. Tanaka for helpfulcomments on the version of this paper that was posted as arXiv:1406.6691v1,the two anonymous referees for useful reports. MRK and JCF acknowledgesupport from NSF CAREER grant AST-0955300. MRK thanks the AspenCenter for Physics, which is supported by NSF Grant PHY-1066293, forproviding the forum during which this work was begun.
References
Adams, F.C., Hollenbach, D., Laughlin, G., Gorti, U., 2004. Photoevap-oration of Circumstellar Disks Due to External Far-Ultraviolet Radia-tion in Stellar Aggregates. ApJ 611, 360–379. doi: , arXiv:arXiv:astro-ph/0404383 .Anderson, D.G., 1965. Iterative procedures for nonlinear integral equations.J. ACM 12, 547–560. URL: http://doi.acm.org/10.1145/321296.321305 , doi: .Benz, W., Ida, S., Alibert, Y., Lin, D.N.C., Mordasini, C., 2014. PlanetPopulation Synthesis, in: Beuther, H., Klessen, R.S., Dullemond, C.P.,Henning, T. (Eds.), Protostars and Planets VI. arXiv:1402.7086 . in press,arXiv:1402.7086.Bertin, G., Lodato, G., 1999. A class of self-gravitating accretion disks. A&A350, 694–704. arXiv:arXiv:astro-ph/9908095 .Bhattacharjee, P., Chaudhury, S., Kundu, S., 2014. Rotation Curve of theMilky Way out to ˜200 kpc. ApJ 785, 63. doi: , arXiv:1310.2659 . 49randenburg, A., Dobler, W., 2002. Hydromagnetic turbulence in computersimulations. Computer Physics Communications 147, 471–475. doi: , arXiv:arXiv:astro-ph/0111569 .Burden, R.L., Faires, J.D., 2011. Numerical Analysis, 9th edition. ThomsonBrooks/Cole.Calef, M.T., Fichtl, E.D., Warsa, J.S., Berndt, M., Carlson, N.N.,2013. Nonlinear Krylov acceleration applied to a discrete ordinatesformulation of the k-eigenvalue problem. J. Comp. Phys. 238, 188–209. URL: , doi: doi:10.1016/j.jcp.2012.12.024 .Cambier, H.J., Smith, D.M., 2013. Prototyping Non-equilibrium Viscous-timescale Accretion Theory Using LMC X-3. ApJ 767, 46. doi: , arXiv:1303.6218 .Colella, P., Woodward, P.R., 1984. The Piecewise Parabolic Method (PPM)for Gas-Dynamical Simulations. Journal of Computational Physics 54,174–201. doi: .Crank, J., Nicolson, P., 1996. A practical method for numerical evaluationof solutions of partial differential equations of the heat-conduction type.Advances in Computational Mathematics 6, 207–226. URL: http://dx.doi.org/10.1007/BF02127704 , doi: .Daniel, J.W., Gragg, W.B., Kaufman, L., Stewart, G.W., 1976. Re-orthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization. Math. Comput. 30, 772–795. URL: ,doi: http://dx.doi.org/10.1090/S0025-5718-1976-0431641-8 .Elmegreen, B.G., 2011. Gravitational Instabilities in Two-component GalaxyDisks with Gas Dissipation. ApJ 737, 10. doi: , arXiv:1106.1580 .Fang, H.r., Saad, Y., 2009. Two classes of multisecant methods for nonlinearacceleration. Numerical Linear Algebra with Applications 16, 197–221.URL: http://dx.doi.org/10.1002/nla.617 , doi: .50letcher, C.A.J., 1991. Computational Techniques for Fluid Dynamics, vol. 1,2nd ed. Springer.Forbes, J., Krumholz, M., Burkert, A., 2012. Evolving Gravitationally Un-stable Disks over Cosmic Time: Implications for Thick Disk Formation.ApJ 754, 48. doi: , arXiv:1112.1410 .Forbes, J.C., Krumholz, M.R., Burkert, A., Dekel, A., 2014. Balanceamong gravitational instability, star formation and accretion determinesthe structure and evolution of disc galaxies. MNRAS 438, 1552–1576.doi: .Gans, P., Gill, J.B., 1984. Smoothing and differentiation of spectro-scopic curves using spline functions. Applied Spectroscopy 38, 370–376. URL: , doi: doi:10.1366/0003702844555511 .Gorti, U., Dullemond, C.P., Hollenbach, D., 2009. Time Evolution of ViscousCircumstellar Disks due to Photoevaporation by Far-Ultraviolet, Extreme-Ultraviolet, and X-ray Radiation from the Central Star. ApJ 705, 1237–1251. doi: , arXiv:0909.1836 .Gorti, U., Hollenbach, D., 2009. Photoevaporation of Circumstellar DisksBy Far-Ultraviolet, Extreme-Ultraviolet and X-Ray Radiation from theCentral Star. ApJ 690, 1539–1552. doi: , arXiv:0809.1494 .Horn, B., Lyra, W., Mac Low, M.M., S´andor, Z., 2012. Orbital Migration ofInteracting Low-mass Planets in Evolutionary Radiative Turbulent Mod-els. ApJ 750, 34. doi: , arXiv:1202.1868 .Hueso, R., Guillot, T., 2005. Evolution of protoplanetary disks: constraintsfrom DM Tauri and GM Aurigae. A&A 442, 703–725. doi: , arXiv:astro-ph/0506496 .Knoll, D.A., Keyes, D.E., 2004. Jacobian-free newton–krylov methods:a survey of approaches and applications. J. Comp. Phys. 193, 357– 397. URL: , doi: http://dx.doi.org/10.1016/j.jcp.2003.08.010 . 51ruijssen, J.M.D., Dale, J.E., Longmore, S.N., 2015. The dynamical evo-lution of molecular clouds near the Galactic Centre - I. Orbital structureand evolutionary timeline. MNRAS 447, 1059–1079. doi: , arXiv:1412.0664 .Krumholz, M.R., Burkert, A., 2010. On the Dynamics and Evolution ofGravitational Instability-dominated Disks. ApJ 724, 895–907. doi: , arXiv:1003.4513 .Liu, B.F., Mineshige, S., Meyer, F., Meyer-Hofmeister, E., Kawaguchi, T.,2002. Two-Temperature Coronal Flow above a Thin Disk. ApJ 575, 117–126. doi: , arXiv:astro-ph/0204174 .Lynden-Bell, D., Pringle, J.E., 1974. The evolution of viscous discs and theorigin of the nebular variables. MNRAS 168, 603–637.Lyra, W., Paardekooper, S.J., Mac Low, M.M., 2010. Orbital Migrationof Low-mass Planets in Evolutionary Radiative Models: Avoiding Catas-trophic Infall. ApJ 715, L68–L73. doi: , arXiv:1003.0925 .Mayer, M., Pringle, J.E., 2007. Time-dependent models of two-phase ac-cretion discs around black holes. MNRAS 376, 435–456. doi: , arXiv:astro-ph/0612751 .McKinney, J.C., Gammie, C.F., 2002. Numerical Models of Viscous Accre-tion Flows near Black Holes. ApJ 573, 728–737. doi: , arXiv:astro-ph/0204045 .McKinney, J.C., Gammie, C.F., 2013. VHD: Viscous pseudo-Newtonianaccretion. Astrophysics Source Code Library. arXiv:1306.015 .Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O.,Zanni, C., Ferrari, A., 2007. PLUTO: A Numerical Code for Com-putational Astrophysics. ApJS 170, 228–242. doi: , arXiv:arXiv:astro-ph/0701854 .Paczy´nsky, B., Wiita, P.J., 1980. Thick accretion disks and supercriticalluminosities. A&A 88, 23–31. 52apaloizou, J.C.B., Stanley, G.Q.G., 1986. The structure and stability ofthe accretion disc boundary layer. MNRAS 220, 593–610.Popham, R., Narayan, R., 1992. Supersonic infall and causality in accretiondisk boundary layers. ApJ 394, 255–260. doi: .Popham, R., Narayan, R., Hartmann, L., Kenyon, S., 1993. Boundary Layersin Pre-Main-Sequence Accretion Disks. ApJ 415, L127. doi: .Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1992. Nu-merical Recipes in C, 2nd Edition. Cambridge University Press, New York.Pringle, J.E., 1981. Accretion discs in astrophysics. ARA&A 19, 137–162.doi: .Reichel, L., Gragg, W.B., 1990. Algorithm 686: Fortran subroutines forupdating the qr decomposition. ACM Trans. Math. Softw. 16, 369–377.URL: http://doi.acm.org/10.1145/98267.98291 , doi: .Saad, Y., Schultz, M., 1986. Gmres: A generalized mini-mal residual algorithm for solving nonsymmetric linear systems.SIAM Journal on Scientific and Statistical Computing 7, 856–869. URL: http://dx.doi.org/10.1137/0907058 , doi: , arXiv:http://dx.doi.org/10.1137/0907058 .Shakura, N.I., Sunyaev, R.A., 1973. Black holes in binary systems. Observa-tional appearance. A&A 24, 337–355.Shu, F.H., 1992. Physics of Astrophysics, Vol. II. University Science Books.Stone, J.M., Mihalas, D., Norman, M.L., 1992. ZEUS-2D: A radiation mag-netohydrodynamics code for astrophysical flows in two space dimensions.III - The radiation hydrodynamic algorithms and tests. ApJS 80, 819–845.Stone, J.M., Norman, M.L., 1992a. ZEUS-2D: A radiation magnetohydro-dynamics code for astrophysical flows in two space dimensions. I - Thehydrodynamic algorithms and tests. ApJS 80, 753–790.53tone, J.M., Norman, M.L., 1992b. ZEUS-2D: A Radiation Magnetohydro-dynamics Code for Astrophysical Flows in Two Space Dimensions. II. TheMagnetohydrodynamic Algorithms and Tests. ApJS 80, 791.Tanaka, T., 2011. Exact time-dependent solutions for the thin accretion discequation: boundary conditions at finite radius. MNRAS 410, 1007–1017.doi: , arXiv:1007.4474 .Toomre, A., 1964. On the gravitational stability of a disk of stars. ApJ 139,1217–1238.Visser, R., Dullemond, C.P., 2010. Sub-Keplerian accretion onto cir-cumstellar disks. A&A 519, A28. doi: , arXiv:1005.1261 .Walker, H., Ni, P., 2011. Anderson acceleration for fixed-point iterations. SIAM Journal on Numerical Anal-ysis 49, 1715–1735. URL: http://epubs.siam.org/doi/abs/10.1137/10078356X , doi: , arXiv:http://epubs.siam.org/doi/pdf/10.1137/10078356X .Willert, J., Taitano, W.T., Kroll, D., 2014. Leveraging Anderson Accelera-tion for improved convergence of iterative solutions to transport systems.J. Comp. Phys. 273, 278–286. Appendix A. Tabulated Rotation Curves
VADER allows arbitrary rotation curves v φ ( r ), and these can be specifiedeither via user-supplied analytic functions, or in the form of a table of ( r, v φ )values. In the latter case VADER generates a rotation curve v φ on the grid viainterpolation, and the potential ψ by integrating the interpolating function.However, the interpolation procedure requires special care to ensure smooth-ness. Equations (1) and (2) involve the second derivative of the torque T ,which is proportional to the rotation curve index β = ∂ ln v φ /∂ ln r . Thusderivatives up to ∂ ln v φ /∂ ln r appear in the evolution equations, and it istherefore highly desirable, for both computational stability and to avoid im-posing artifacts on the solution, that the interpolating function constructedto approximate a table of ( r, v φ ) values have at least three continuous deriva-tives. 54o achieve this aim, VADER constructs tabulated rotation curves usingbasis splines (b-splines). In the b-spline method, the domain of the functionto be fit is broken up into a set of intervals, separated by breakpoints, andone must choose the degree D of the fit, the number of breakpoints B , andtheir locations. For the choice of locations, we use the method of Gans andGill (1984), who show that errors in the resulting fit are minimized if thebreakpoints are distributed evenly in the size of the the interval weighted bythe square root of the function being fit. Let ( r n , v φ,n ) be our table of N input data, ordered from n = 0 . . . N − r n increasing monotonically,and let b m be the location of the m th breakpoint, m = 0 . . . B −
1. We set b = r and b B − = r N − , and we assign the remaining breakpoints b m viathe following algorithm. Let S T ≡ B + 1 N − (cid:88) n =0 v / φ,n dx n , (A.1)where dx n = (1 /
2) log( r n +1 /r n − ) for n (cid:54) = 0 , N − dx = log( r /r ), and dx N − = log( r N − /r N − ). Starting from a breakpoint b m located at datavalue r n m , we set the position of the next breakpoint to b m +1 = r n m +1 , where n m +1 is the smallest index for which n m +1 (cid:88) n = n m +1 v / φ,n dx n ≥ S T . (A.2)This assures that the breakpoints are as uniformly distributed as possiblefollowing the criterion of Gans and Gill (1984). Once the breakpoint loca-tions are chosen, the function and its derivatives can be constructed by thestandard basis spline method.To demonstrate and evaluate the performance of this capability, we per-form two tests. For the first, we take the Paczy´nsky and Wiita (1980) ap-proximation for a black hole potential ψ PW = GMr − r g , (A.3) If the number of points N is very small, it is conceivable that, due to the discreteplacement of the data points r n , this algorithm might result in not all breakpoints b m being assigned before we reach the end of the array. In this case we end up with the lastbreakpoint interval being of size 0, i.e., b B − = b B − . However, this presents no obstaclesto the remainder of the algorithm. r g = 2 GM/c is the horizon radius, and generate a table of rotation ve-locities by analytically evaluating v φ = dψ PW /dr at points from 1 . r g − . r g ,spaced in units of 0 . rg g . We then use this table to generate a 6th-order b-spline fit using 15 breakpoints, and from that fit compute the potential ψ and the rotation curve index β on a grid of 512 points logarithmically spacedfrom r = 2 − r g . Figure A.12 shows the results of the fit as compared withthe analytic values for ψ , v φ , and β . Clearly the b-spline reconstruction isexcellent.For the second test we use a much noisier data set: a compilation of dataon the rotation curve of the Milky Way from Bhattacharjee et al. (2014),using their model where the Sun is 8.5 kpc from the Galactic center andthe rotation velocity at the Solar circle is 220 km s − . Figure A.13 showsthe data and a variety of VADER b-spline reconstructions of v φ , ψ , and β on a logarithmic grid of 512 points uniformly spaced from 0 . −
180 kpc.In this case it is clear that the rotation curve and potential are very wellreconstructed and that fits to them do not depend strongly on the choiceof degree D and number of breakpoints B , except that D ≥ β does depend at least somewhaton these choices. For this particular data set, there appears to be no valueof D that both guarantees that the derivative of β be continuous and that β itself remains within physically-reasonable values (roughly − . β , prune the data set more carefully, or use a more sophisticated fittingprocedure. 56 v φ , ψ , β Analytic v φ ψβ Fit v φ ψβ r/r g -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 | E rr o r | Error > 0Error < 0
Figure A.12: B-spline reconstruction of the rotation curve for a Paczy´nsky and Wiita(1980) potential using a fit of degree D = 6 with B = 15 breakpoints. In the upper panel,we show the exact analytic values of v φ , ψ , and β (solid lines) and the numerical fits (datapoints, only every 16th point shown for clarity). All quantities are plotted in units where GM = r g = 1, and the gauge of the potential set so that ψ → r → ∞ . In thebottom panel, we plot the relative error of the b-spline reconstruction versus r , defined asError = ( v φ, fit − v φ, exact ) /v φ, exact , and similarly for ψ and β . v φ [ k m s − ] Data D =2 D =4 D =3 , B =15 D =3 , B =8 ψ / ( k m s − ) -1 r [kpc]32101 β Figure A.13: B-spline reconstruction of the rotation curve of the Milky Way. The toppanel shows data from Bhattacharjee et al. (2014) (black points with error bars), togetherwith b-spline fits of degree D = 2, 3, and 4 (blue lines). The D = 2 and D = 4 fitsuse 15 breakpoints each, while for D = 3 we show models with both B = 8 and B = 15breakpoints. The middle panel shows the potential, with a gauge chosen so that ψ = 0 atthe edge of the grid. The bottom panel shows β ..