Issue 
A&A
Volume 545, September 2012



Article Number  A152  
Number of page(s)  16  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201219557  
Published online  25 September 2012 
A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO^{⋆} code
^{1} Dipartimento di Fisica GeneraleUniversitá di Torino, via Pietro Giuria 1, 10125 Torino, Italy
email: mignone@ph.unito.it
^{2} CEA Irfu, SAP, Centre de Saclay, 91191 GifsurYvette, France
^{3} Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
email: flock@mpia.de
^{4} Institute for Astronomy and Astrophysics, Section Computational Physics, Eberhard Karls Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
email: matthias.stute@tat.physik.unituebingen.de; stefan.kolb@tat.physik.unituebingen.de
^{5} Consorzio Interuniversitario CINECA, via Magnanelli 6/3, 40033 Casalecchio di Reno, Bologna, Italy
Received: 8 May 2012
Accepted: 29 June 2016
Context. Explicit numerical computations of hypersonic or superfast differentially rotating disks are subject to the timestep constraint imposed by the Courant condition, according to which waves cannot travel more than a fraction of a cell during a single timestep update. When the bulk orbital velocity largely exceeds any other wave speed (e.g., sound or Alfvén), as computed in the rest frame, the time step is considerably reduced and an unusually large number of steps may be necessary to complete the computation.
Aims. We present a robust numerical scheme to overcome the Courant limitation by improving and extending the algorithm previously known as FARGO (fast advection in rotating gaseous objects) to the equations of magnetohydrodynamics (MHD) using a more general formalism. The proposed scheme conserves total angular momentum and energy to machine precision and works in cartesian, cylindrical, or spherical coordinates. The algorithm has been implemented in the next release of the PLUTO code for astrophysical gasdynamics and is suitable for local or global simulations of accretion or protoplanetary disk models.
Methods. By decomposing the total velocity into an average azimuthal contribution and a residual term, the algorithm approaches the solution of the MHD equations through two separate steps corresponding to a linear transport operator in the direction of orbital motion and a standard nonlinear solver applied to the MHD equations written in terms of the residual velocity. Since the former step is not subject to any stability restriction, the Courant condition is computed only in terms of the residual velocity, leading to substantially larger time steps. The magnetic field is advanced in time using the constrained transport method in order to fulfill the divergencefree condition. Furthermore, conservation of total energy and angular momentum is enforced at the discrete level by properly expressing the source terms in terms of upwind Godunov fluxes available during the standard solver.
Results. Our results show that applications of the proposed orbitaladvection scheme to problems of astrophysical relevance provides, at reduced numerical cost, equally accurate and less dissipative results than standard timemarching schemes.
Key words: methods: numerical / accretion, accretion disks / protoplanetary disks / magnetohydrodynamics (MHD) / turbulence
Freely distributed at http://plutocode.ph.unito.it.
© ESO, 2012
1. Introduction
The physics of accretion flows has received a great deal of attention over the past few decades, and several numerical investigations have contributed to broaden our current understanding of the problem in its diverse aspects. Representative astrophysical scenarios involve magnetically driven turbulence in accretion disks around compact objects, which are typically modeled via a local approach (the socalled shearingbox model of Hawley et al. 1995) or global disk simulations, that have recently become more appealing owing to increases in the available computational power (see, for instance Beckwith et al. 2011; Flock et al. 2011; Sorathia et al. 2012). Likewise, in the context of planet formation, simulations of protoplanetary disks have become an important tool for the study of the dynamics of the gas and its impact on either particles or planets (see Kley et al. 2009; Uribe et al. 2011, and references therein).
A common ingredient of all models is a highly supersonic or superfast, rotating, sheared flow. In standard explicit numerical calculations, the time step is restricted, for stability reasons, by the CourantFriedrichsLewy (CFL) condition (Courant et al. 1928), which is formally stated as (1)where Δl is the cell length, λ_{max} is the fastest signal speed, and C_{a} – the Courant number – is a limiting factor typically less than one. In essence, the CFL condition given in Eq. (1) prevents any wave from traveling more than a fraction of a grid cell. This leads to a drastic reduction in the time step, whenever the orbital speed becomes considerably greater than either the sound or fast magnetosonic velocities, and results in excessively long computations. A typical occurrence is encountered in accretion or protoplanetary disks where the Keplerian velocity is the dominating flow speed.
In the context of hydrodynamics, Masset (2000) presented the first fast Euleriantransport algorithm for differentially rotating disks (FARGO) in which the timestep is limited by the velocity fluctuations around the mean orbital motion rather than by the total azimuthal velocity itself. The algorithm is based on the simple recognition that if the orbital velocity is written as the sum of a constant average term plus a residual, the evolution operator can be implemented as the composition of a linear transport operator and a nonlinear step involving the residual velocity. Since the former step yields a uniform shift along the orbital direction, the CFL condition is determined during the nonlinear step by the magnitude of the residual velocity. The algorithm may be viewed as a temporary change to the local corotating frames placed at different radial positions of the differentially rotating disk. The advantage of using such a strategy is twofold: on the one hand, it reduces numerical dissipation as the Keplerian flow is now treated as a mean flow in the equations and, on the other hand, it allows larger time steps to be taken thus reducing the computational costs (as long as the residual velocity remains subsonic).
The ideas and concepts behind FARGO were extended to magnetohydrodynamics (MHD) by Johnson et al. (2008b) and Stone & Gardiner (2010) in the context of shearingbox model which is a local approximation describing a small disk patch embodied by a cartesian box corotating with the disk at some fiducial radius. Because of the difficulties inherent in simulating an entire disk, local shearingbox models have largely contributed to much of what is presently known about the magnetorotational instability and its implication for the transport of angular momentum in disks. Nevertheless, several authors have stressed, over the past decade, the importance of global disk models, specially in the context of magnetorotational instability, e.g., Armitage (1998), Hawley (2000), Fromang & Nelson (2006) and Regev & Umurhan (2008). Owing to the increased numerical power, highresolution simulations of global disk models are now becoming amenable to investigation. In this context, a first implementation of the FARGO scheme for MHD was presented by Sorathia et al. (2012) as part of the Athena code (Stone et al. 2008) in cylindrical geometry.
Here we propose, using a somewhat different and more general formalism, an extension of the original FARGO scheme to the equations of MHD for different systems of coordinates in the context of Godunovtype schemes. The algorithm differs from the original one in that it uses a fully unsplit integration scheme when solving the MHD equations in the residual velocity. Source terms arising from the velocity decomposition are treated carefully in order to restore the conservation of total energy and total angular momentum. The magnetic field is evolved via the constrainedtransport (CT) scheme, which preserves the divergencefree condition to machine accuracy at all times. The proposed FARGOMHD scheme is implemented in the PLUTO code for astrophysical gasdynamics (Mignone et al. 2007, 2012) and will be available with the next code release (4.0). The algorithm has been fully parallelized so that domain decomposition can be performed in all three coordinate directions, thus allowing an efficient use of a large number of processors.
The paper is organized as follows. In Sect. 2, we briefly review the basic MHD equations and the timestep limitations imposed by standard explicit timestepping methods (Sect. 2.1). The new FARGOMHD scheme is presented in Sect. 2.2 and consists of i) a standard nonlinear step (Sect. 2.3) formulated so as to preserves total energy and total angularmomentum conservation (Sect. 2.3) and ii) a linear transport step (Sect. 2.4) where fluid quantities are simply advected in the direction of orbital motion. Numerical benchmarks and astrophysical applications are presented in Sect. 3, while conclusions are drawn in Sect. 4.
2. Description of the FARGOMHD scheme
2.1. Standard approach and limitations
We begin by considering the ideal MHD equations written as a nonlinear system of conservation laws where ρ is the mass density, v is the fluid velocity, Φ is the gravitational potential, B is the magnetic field, p_{t} = p + B^{2}/2 is the total pressure accounting for thermal (p), and magnetic (B^{2}/2) contributions. The total energy density E is given by (6)where an ideal equation of state (EoS) with specific heat ratio Γ has been used. Dissipative effects are neglected for the sake of simplicity, although they can be easily incorporated in this framework.
In the usual finitevolume approach, Eq. (2) to Eq. (5) are discretized on a computational mesh spanned by the grid indices i, j, and k corresponding to the location of a given cell in a threedimensional (3D) coordinate system. We label the generic unit vector along a given axis with and consider, in what follows, cartesian (), cylindrical polar (), and spherical () coordinates. Individual cells have spatial extents given by Δl_{d,ijk}, where d labels the direction.
The time step Δt is determined by the CFL condition in Eq. (1), the precise form of which can depend on the employed timestepping scheme. In the PLUTO code, one can take advantage of either the cornertransport upwind (CTU) method of Colella (1990) and Gardiner & Stone (2008) or resort to a RungeKutta (RK) discretization, where the spatial gradients are treated as the righthand side of an ordinary differential equation (method of lines). Both algorithms are dimensionally unsplit but stable under somewhat different Courant conditions. Here we consider the 6solve CTU and the secondorder RungeKutta (RK2) scheme, which have the same number of Riemann problems per cell per step. If N_{dim} = 2,3 is the number of spatial dimensions, one has (Beckers 1992; Toro 1999) (7)where v_{d} and c_{f,d} are the fluid velocity and fast magnetosonic speed in the direction d. For highly superfast, gridaligned flows, Eq. (7) shows that both CTU and RK2 yield comparable time steps in twodimensions, while, in threedimensions, RK2 is expected to give time steps that are approximately twice as large.
Nevertheless, numerical computations of advectiondominated flows for which v_{d} ≫ c_{f} may result in small time steps yielding unusually long computations and eventually suffering from an excess of numerical dissipation. This kind of scenario worsens in the case of a Keplerian flow in polar or spherical coordinates, since the cell length becomes increasingly smaller towards the inner radius where the orbital flow velocity is faster.
2.2. An orbital advection scheme for MHD
To overcome the limitations outlined in the previous section, we now assume that the fluid velocity may be decomposed as v = v′ + w, where w is a solenoidal velocity field and v′ is the fluctuation or residual. Equations (2) through (5) may then be written as where I is the identity matrix and E′ is the residual energy density, defined as (12)and the two additional source terms in the momentum and energy equations may be written, after some algebra, as (13)and (14)Explicit expressions for these terms in different systems of coordinates can be found in Appendix A.
We note that, Eqs. (8) through (11) have a general validity since no assumption has been made about the magnitude of the residual velocity v′. However, neither the residual momentum nor the residual energy density are conserved quantities because of the appearance of the additional source terms and given by Eqs. (13) and (14). Since nothing has of course changed at the continuous level, both the total energy and momentum must still be conserved. The situation is different, however, at the discrete level where differential vector identities are satisfied only within the truncation error of the scheme. This is discussed in more detail in Sect. 2.3.1
If we look at the individual components of the system of Eqs. (8) through (11), each evolution equation has the form (15)where q ∈ (ρ,ρv′,E′,B), F_{q} is the usual MHD flux written in terms of the residual velocity v′. The source term S_{q} appearing on the righthand side of Eq. (15) accounts for several contributions that include gravity, the additional term (or ), and, in the case of curvilinear coordinate systems, geometrical factors arising from differential operators. The last term on the lefthand side of Eq. (15) describes the linear transport of q with advection speed w. If w coincides with the average azimuthal velocity, this term has the simple effect of pushing fluid elements along their orbit.
An effective approach to solve Eq. (15) is to use operator splitting and separate the solution into a standard nonlinear step that does not include the linear advection term w·∇q(16)and a linear transport step corresponding to the solution of (17)The two operators are separately described in Sects. 2.3 and 2.4, respectively.
When the orbital motion coincides with one of the coordinate axes, e.g. , the exact solution of (17) can be formally expressed as q(y,t^{n + 1}) = q(y − wΔt^{n},t^{n}), which corresponds to a noninteger translation of a row of computational zones. Since this operation is decomposed into an integer shift plus a remainder, the linear transport step is unconstrained and can be carried out for arbitrarily large time steps.
The time step is thus ultimately determined by the Courant condition during the nonlinear step, that is, by solving the regular MHD equations (plus additional source terms) written in terms of the residual velocity v′. Provided that , the CFL condition in Eq. (7) now results in appreciably larger time steps, since the dominant background orbital contribution has now been removed from the computation.
2.3. Standard nonlinear step
During the standard nonlinear step, we solve Eqs. (8) through (11) without the linear advection operator (w·∇). The resulting system is equivalent to the regular MHD equations written in terms of the residual velocity v′ = v − wplus the additional source terms and given, respectively, by Eqs. (13) and (14).
This entitles us to employ the formalism already developed for the solution of the MHD equations using any stable conservative finitevolume or finitedifference discretization available with the PLUTO code. Accordingly, we update each conserved quantity q using the general building block (18)which is the discrete version of Eq. (16). In Eq. (18), the flux F_{q} follows from the solution of Riemann problems at zone interfaces, while the divergence term is expressed as the sum of twopoint difference operators in each direction d. Dropping the index q for simplicity and considering a generic flux F with components , we use (19)where F_{d, + } and F_{d, − } are the fluxes through the upper (+) and lower ( − ) cell boundaries with surface normal , while A_{d} and are the interface areas and cell volume.
The magnetic field is evolved using the constrained transport formalism (see, e.g., the papers by Mignone et al. 2007; Flock et al. 2010).
2.3.1. Conservative formulation
As anticipated in Sect. 2.2, the source term S_{q} contains both pointwise contributions (e.g. gravitational or curvilinear terms) and terms involving the derivatives of the orbital velocity w. Pointwise source terms not involving spatial derivatives are added to the righthand side of the equations in a standard explicit way. Conversely, source terms contributing to both the azimuthal momentum component and the energy equation require some additional considerations. To this end, we consider the net change in the total azimuthal momentum during a single timestep update. Using, for simplicity, cartesian coordinates and neglecting gravity, a straightforward combination of the equations of continuity and the ycomponent of momentum yields (20)Although the sum of the last two terms on the righthand side is equal to ∇·(wF_{ρ}) at the continuous level, this may not necessarily hold in a discrete sense. At the numerical level one should indeed expect  w∇·F_{ρ} + ρv′·∇w − ∇·(wF_{ρ})  = O(Δx^{s}), where O(Δx^{s}) is the truncation error of the scheme. As a consequence, total (linear or angular) momentum and, similarly, total energy density will be conserved only at the truncation level of the scheme.
We, instead, seek a numerical discretization that allows us to restore the exact conservation of total linear (for cartesian geometry) or angular (for polar grids) momentum and energy, also at the numerical level. For this purpose, we return to the derivations given in Eqs. (2) − (5) to (8) − (11) and note that the source term, as given in Eq. (13), of the momentum equation simply results from the algebraic manipulation of (21)where the last term is identically zero in cartesian coordinates and only contributes to the radialmomentum source term in polar and spherical coordinates. Similarly, one can show, after some algebra, that the source term in Eq. (14) of the energy equation results from (22)This suggests that Eqs. (21) and (22) may be conveniently expressed in terms of the density and ymomentum upwind fluxes already at disposal during the conservative update. For instance, we update density, (residual) azimuthalmomentum, and (residual) energy as (23)(24)(25)which shows that, by adding the product of Eq. (23) and w to Eq. (24), one obtains the conservative update of the total momentum with corresponding flux F_{my} + wF_{ρ}. Since similar arguments hold by combining the density and momentum equations with Eq. (25), the resulting discretization enforces the conservations of both total momentum and total energy in both a continuous and a numerical sense. The corresponding expressions for polar and spherical coordinates are a straightforward extension of these equations leading to both total angular momentum and total energy conservation. We report them in Sects. A.2 and A.3, respectively. We point out that the importance of a conservative formulation of the equations has already been recognized by other authors, e.g. Kley (1998), who has shown that the inclusion of nonintertial forces as a source term may lead to erroneous results.
2.4. Linear transport step
As anticipated, the solution of the linear transport equation in Eq. (17) consists of a uniform shift  wΔt  of the conserved variable profiles in the direction of orbital motion. If the flow is predominantly aligned with one of the coordinate axes, say , we can write the solution of Eq. (17) for t ∈ [t^{n},t^{n + 1}] as q(y,t) = q(y − w(t − t^{n}),t^{n}). Extensions to polar and spherical coordinates are straightforward, provided that y → φ and wt → Ωt.
Integration of Eq. (17) for a time step Δt and over a cell size Δy gives (26)where q^{n}(ξ) = q(ξ,t^{n}) and w does not depend on y. The integrals on the righthand side of Eq. (26) can be converted into a finite sum corresponding to an integer shift of cells m = floor(wΔt/Δy + 1/2) plus a fractional remainder δy = wΔt − mΔy. The final result may be cast as a twopoint flux difference scheme (27)where j_{m} = j − m and (28)is the upwind numerical flux. We remark that the finite summations corresponding to an integer cell shift (implicitly contained in the integral of Eq. (26)) do not need to be explicitly computed since most terms cancel out when taking their difference. We also note that, since  δy  < Δy by construction, the scheme is always stable regardless of the choice of Δt.
The fluxes given by Eq. (28) may be computed to an arbitrary order of accuracy by assuming a piecewise polynomial representation of the data inside the cell. If, for instance, we assume a piecewise linear distribution of q inside the zone, then the integral in Eq. (28) takes the form (29)where Δq_{j} is computed using a standard slope limiter. Equation (29) corresponds to the wellknown secondorder MUSCLHancock scheme of van Leer (1984). Similarly, using a piecewise parabolic distribution inside the cell we have (30)where q_{j, ± } are third (or higher) order rightmost and leftmost interpolated values with respect to the cell center, δq_{j} = q_{j, + } − q_{j, − } and δ^{2}q_{j} = q_{j, + } − 2q_{j} + q_{j, − }. This is essentially the PPM scheme for a linear advection equation (Colella & Woodward 1984) and is our default choice, unless stated otherwise.
2.4.1. Magnetic field transport
In the constrained transport formalism, the three components of magnetic field are discretized on a staggered mesh and evolved as surface averages placed at the different zone faces to which they are orthogonal. In this sense, we locate the components of B by means of a halfinteger subscript, that is, , and to denote the magnetic field components centered on the x, y, and z faces, respectively.
The evolution is carried out using a discrete version of Stokes’ theorem applied to Eq. (11), where a line integral of the electric field is properly evaluated at zone edges. This guarantees that the divergencefree condition of the magnetic field is also maintained to machine precision during the linear transport step. During the transport step (i.e. when v′ = 0), this leads to the following discretization of Eq. (11) where ℰ^{x} = −∫wB_{z} dt and ℰ^{z} = ∫wB_{x} dt are the x and z components of the timeintegrated electromotive force computed, in analogy with the results of the previous section, as (34)with computed from Eq. (28) using q = B_{x} and (35)with computed from Eq. (28) with q = −B_{z}. In the previous equations, we have dropped, for simplicity, the full subscript notations when redundant and kept only the halfincrement notation to represent the different magnetic and electric field components.
We point out that the updating formulas for B_{x} and B_{z} could be implemented in exactly the same way as done for Eq. (27) since the differences in the fluxes along the ydirection, defined by Eqs. (31) and (33), lead to the cancellation of most terms leaving only the upstream values. On the other hand, cancellation does not occur when updating B_{y}, since Eq. (32) involves differences of terms along the x and z directions corresponding to the winding of the field under the action of the shear. For this reason, the full summation must now be retained in Eqs. (34) and (35) when evaluating the timeintegrated electromotive force.
For the sake of completeness, we also report the corresponding update expressions for other geometries. In polar coordinates (R,φ,z), we have where ℰ^{R} = −^{∫}ΩB_{z} dt, ℰ^{z} = ^{∫}ΩB_{R} dt (with Ω = w/R) are computed similarly to Eqs. (34) and (35).
Likewise, in spherical coordinates (r,θ,φ) we have where ℰ^{r} = ∫ΩB_{θ} dt, ℰ^{θ} = −∫ΩB_{r} dt, and Ω = w/(rsinθ) are computed similarly to Eqs. (34) and (35).
2.5. Expected speedup
In general, the expected timestep gain is problemdependent but can also be affected by several factors such as the magnitude of the orbital velocity, the grid geometry, and the cell aspect ratio. An estimate of the speedup may be directly computed from Eq. (7) by taking the ratio of the time step obtained with orbital advection to that obtained without. Specializing, for example, to polar coordinates and RK2 timestepping, one can estimate a timestep increase (42)If the characteristic signal fluctuations are of comparable magnitude, and are such that λ′ ≪ w, the previous expression simplifies to (43)where M = w/λ′ should be of the same order as the fast magnetosonic Mach number. The previous estimate shows that a sensible choice of the grid resolution has a direct impact on the expected gain: a finer resolution in the azimuthal direction (RΔφ ≪ ΔR) favors a larger gain, whereas cells with a smaller radial extent may become less advantageous.
2.6. Parallel implementation
In parallel computations, the simulation domain is decomposed into smaller subgrids, which are then solved simultaneously by several processing units. For timeexplicit calculations, only boundary data stored in the ghost zones have to be exchanged between neighboring processors. With FARGOMHD, however, this approach presents some difficulties since the linear transport step may involve shifts of an arbitrary number of zones along the orbital direction. This implies that data values could in principle be exchanged between all processors lying on the same row of a parallel domain decomposition, thus involving many more communications than a regular exchange of ghost zones between adjacent processors.
In the tests and applications presented here, we observed that the maximum zone shift m_{max} hardly ever exceeds the typical grid size N_{φ} of a single processor. As efficient parallel applications require N_{φ} ≳ 8 − 16, we can safely assume that the condition m_{max} ≤ N_{φ} is virtually always respected and does not represent a stringent requirement for most applications. In such a way, parallel communications are performed as a cyclic shift between neighboring processes only along the orbital direction and the computational overhead becomes approximately 1 + m_{max}/n_{g} times larger than a regular boundary call (where n_{g} is the number of ghost zones) in any direction. Notice also that, when the orbital velocity does not change sign across the domain, the amount of subtask communication can be halved since information always travels in the same direction and data values need to be transferred from one processor to the next following the same pattern.
3. Numerical benchmarks
We present a number of hydrodynamical and magnetohydrodynamical test problems where the proposed orbital advection algorithm is directly compared with the standard traditional integration method. Unless stated differently, we make use of the ideal EoS with Γ = 5/3 and the PPM method, given by Eq. (30), is the default interpolation during the linear transport step.
3.1. Vortex dynamics in a twodimensional Keplerian flow
We begin by considering the dynamical evolution of an anticyclonic vortex in a twodimensional (2D) Keplerian disk using polar coordinates (R,φ). The initial background state is defined by constant density and pressure equal respectively to ρ = 1 and p = 1/(ΓM^{2}), where M = 10 is the Mach number at R = 1. The disk rotates with angular velocity Ω(R) = R^{ − 3/2} under the influence of a pointmass gravity Φ = −1/R and fills the computational domain defined by 0.4 ≤ R ≤ 2 and 0 ≤ φ ≤ 2π. With these definitions, the disk scaleheight is given by H(R) = c_{s}/Ω(R), where c_{s} = 1/M. A circular vortex is superposed on the mean Keplerian flow and initially defined by (44)where h = H(1)/2 is the size of the vortex in terms of the local scaleheight at R = 1, κ = −1 is the vortex amplitude, and x = Rcosφ − R_{0}cosφ_{0} and y = Rsinφ − R_{0}sinφ_{0} are the cartesian coordinates measured from the center of the vortex initially located at R_{0} = 1 and φ_{0} = π/4.
We solve the hydrodynamical equations (no magnetic field) using the secondorder RungaKutta timestepping with C_{a} = 0.4 using linear reconstruction and follow the system until t = 200π, i.e. for 100 revolutions of the ring at R = 1. Computations are carried out with and without orbital advection at three different resolutions corresponding to 256 × 1024 (low), 512 × 2048 (medium), and 1024 × 4096 (high), yielding approximately square cells in the proximity of the vortex. To enforce conservation, we ensure that no net flux is established across the domain and therefore choose reflective boundary conditions at the inner and outer radial boundaries and impose periodicity along the φ direction.
Fig. 1 Local potential vorticity for the hydrodynamical vortex test after 20 orbits at the resolution of 1024 × 4096. Top panel: standard integration. Bottom panel: orbital advection. 
After a few revolutions, the vortex experiences a nonlinear adjustment to its final persistent structure. This process is accompanied by the emission of spiral density waves, which radiate away the energy excess of the initial state (Bodo et al. 2007). Figure 1 shows the local potential vorticity, or vortensity, defined as after 20 revolutions using the standard integration algorithm (top) and the proposed FARGO scheme (bottom) at the highest resolution. The two snapshots are in good agreement, although the simulation using orbital advection required only ~37 000 steps compared to the ~444 000 steps required by the standard computation. We found indeed that the time step increased, on average, by a factor between ~11 and ~12 with orbital advection, as reported in Table 1.
Average time step for the hydrodynamical vortex test problem at different resolutions using the standard scheme (2nd Col.) and with orbital advection (3rd Col.).
Fig. 2 Time evolution of the vortensity minimum (top, normalized to its initial value) and relative variation in the total integrated vorticity (bottom) for the hydrodynamical vortex problem. Dotted, dashed, and solid lines refer to computations obtained at the resolutions of 256 × 1024, 512 × 2048, and 1024 × 4096 grid zones, respectively. Black and red colors correspond to the the standard integration algorithm and FARGOMHD, respectively. 
As an indicative measure of the dissipation properties of the scheme, we compare, in the top panel of Fig. 2, the time history of the local vortensity minimum for the selected grid sizes. An increase in the grid resolution favors a slower vortex decay leading to the formation of stable, longlived vortex structures (Bodo et al. 2007). The same effect is obtained, at an equivalent resolution, by employing orbital advection.
Similarly, we inspected the relative variation in the total integrated vorticity presented as a function of time in the bottom panel of Fig. 2. Strictly speaking, we note that vorticity in a 2D compressible flow is a conserved quantity only if the fluid is barotropic, i.e., if p = p(ρ). Nevertheless, our results indicate that the amount of generated vorticity decreases from ≈ 18% (at the lowest resolution) to ≈ 10% (at the largest) using the standard integration scheme (black lines), while it remains smaller than 1% when orbital advection is employed (red lines).
Lastly, we verified that both total energy (including the gravitational contribution) and angular momentum are conserved to machine accuracy, as expected from the conservative formulation of our scheme.
3.2. Magnetohydrodynamical vortex in a shear flow
Fig. 3 Evolution of the magnetic pressure distribution, ( × 100), for the MHD vortex test computed using the standard MHD update (top row) and with FARGOMHD (bottom row) using 512^{2} zones. From left to right, we show the evolution from time t = 0 to t = 3. 
In the next example, we investigate the stretching of a 2D MHD vortex in a superfast sheared flow using cartesian coordinates. The initial condition is similar to the one used by Mignone et al. (2010) and consists of a uniformdensity (ρ = 1), background shear flow with velocity profile v_{y} = Mtanh(x)/2, where M is the sonic Mach number. The magnetic field and pressure distributions are, respectively, given by (45)where r^{2} = x^{2} + y^{2} and, for the present example, we adopt μ = 0.02.
We choose the square region x,y ∈ [ − 5,5] as our computational domain and impose periodic boundary conditions in the ydirection while applying reflecting conditions in the xdirections. These choices ensure that no net flux is established across the domain boundaries so that density, momentum and energy are globally conserved. We note that, in the absence of shear (for M = 0), the previous configuration is an exact solution of the ideal MHD equations, representing a circular magnetic fieldloop supported against a pressure gradient. In the present context, however, we choose M = 20 and follow the stretching of the initial loop configuration as time advances. Computations are carried on 512^{2} grid zones using the PPM+CTU scheme with Courant number C_{a} = 0.8.
Results with and without the FARGOMHD algorithm are shown in Fig. 3, where we compare the magnetic pressure distributions for t = 0,1,2,3. The two solutions look indistinguishable, as also visible from a onedimensional cut of the ycomponent of magnetic field at t = 3, shown in the left panel of Fig. 4. A plot of the time step is shown in the right panel of the same figure, where one can see that orbital advection attains a larger gain ( ≳ 9) during the first four revolutions and decreases to ≈ 4.5 towards the end. In terms of CPU time, the standard integration took approximately 2 hours and 31 minutes whereas only 26 min were required using the orbital advection scheme. This gives, on average, a speedup of ≈ 5.8.
To check whether total energy and angular momentum are conserved, we repeated the integration for a much longer time, corresponding to 100 revolutions, and added a random horizontal velocity component of the form (46)where ℛ is a random number between 0 and 1. We compare, in Fig. 5, the conservative and the nonconservative variants of the orbital advection algorithm by plotting the volumeaverage of the y component of momentum and the relative variation in the total energy as time advances. As expected, the former conserves momentum and energy to machine accuracy, while the latter exhibits increasing deviations from zero as the computation proceeds.
Fig. 4 Top panel: horizontal cuts at y = 0 of the y component of magnetic field at t = 3 for the MHD vortex test in cartesian coordinates using the standard integration method (solid line) and FARGOMHD (plus signs). Bottom panel: timestep variation during the computation. Solid and dashed lines correspond to standard and FARGOMHD computations, respectively. 
Fig. 5 Temporal evolution of the volumeaveraged y component of momentum (black) and fractional average energy variation (green) for the MHD vortex with random perturbation using the nonconservative (dashed lines) and conservative (solid line) versions of the orbital advection scheme. 
3.3. Shearingbox models
The shearingbox approximation (Hawley et al. 1995) provides a local cartesian model of a differentially rotating disk. By neglecting curvature terms, one focuses on the evolution of a small rectangular portion of the disk in a frame of reference corotating with the disk at some fiducial radius R_{0} where the orbital frequency is Ω_{0}. In this approximation, the orbital motion is described by a linear velocity shear of the form , where (47)is a local measure of the differential rotation.
The largescale motion of the disk is described by assuming that identical boxes slide relative to the computational domain in the radial direction, a requisite implemented by the shearingbox boundary condition, which enforce sheared periodicity. For any flow quantity q not containing v_{y}, this is formally expressed by (48)where L_{x} is the radial (x) extent of the domain. The implementation of the shearingbox boundary conditions is done similarly to Gressel & Ziegler (2007) and employ the same techniques described in this paper.
We report, for the sake of completeness, the basic equations in Appendix A.1.1 and refer the reader to Hawley et al. (1995), Gressel & Ziegler (2007) and Stone & Gardiner (2010) for a thorough discussion. In the tests below, we employ an isothermal EoS, assume a Keplerian flow (q = 3/2), and adopt the unsplit CTU+PPM scheme with a Courant number C_{a} = 0.45.
3.3.1. Compressible MHD shearing waves
As a first benchmark, we consider the evolution of compressible, magnetohydrodynamical shearingwaves in an isothermal medium using the shearingbox model. These waves can be regarded as the extension of fast and slow modes to a differentially rotating medium and provide the natural basis for wave decomposition in the linear theory of rotating MHD shear flows (Johnson 2007).
We follow the same configuration used by Johnson et al. (2008b) and Stone & Gardiner (2010) and consider a cartesian domain of size [ − 1/4,1/4] ^{3} with N × N/2 × N/2 grid zones. The initial conditions, the details of which may be found more precisely in Johnson et al. (2008a), consists of a background equilibrium state about which a plane wave perturbation is imposed (49)where is the background shear flow, H = c_{s}/Ω_{0} is the scale height, and k = 4π/H( − 2,1,1) is the initial wavenumber, whereas (50)are the perturbation amplitudes of density and vector potential (ϵ = 10^{6}). Magnetic field is initialized from B = ∇ × A. We set the isothermal sound speed to c_{s} = 1, while the orbital frequency is Ω_{0} = 1.
Fig. 6 Azimuthal magneticfield perturbation as a function of time for the compressible MHD shearingwave problem. The exact solution is plotted as a solid black line and the dashed lines correspond to computations obtained at the resolutions of N_{z} = 8 (blue), 16 (green), and 32 (red) zones in the vertical direction. 
In Fig. 6, we compare, at different resolutions N = 8,16,32, the temporal evolution of the azimuthal field perturbation (51)where is the timedependent wavenumber and is the analytical expectation, kindly provided to us by G. Mamatsashvili. Our results are in excellent agreement with those of Johnson et al. (2008b) and Stone & Gardiner (2010) showing that at the highest resolution of 32 zones per wavelength, the solution has converged to the semianalytical prediction.
3.3.2. Nonlinear MRI
Fig. 7 Density distribution in the shearingbox computation after 40 rotation periods. Top and bottom panel refers to the model with and without FARGOMHD, respectively. 
Fig. 8 Evolution of the MRI in the shearingbox model during the first 40 orbits. From top to bottom, the three panels show, respectively, the volumeintegrated magnetic energydensity, the Maxwell stresses, and the time increment used in the computations without (solid line) and with FARGOMHD (dashed line). 
In the next example, we investigate the nonlinear evolution of the magnetorotational instability (MRI) in the shearingbox model. The initial condition consists of a uniform orbital motion in the y direction of constant density ρ = 1 and zero netflux magnetic field in the vertical direction (52)where β = 400 and c_{s} is the isothermal sound speed. Following Gressel & Ziegler (2007), Johnson et al. (2008b), and Stone & Gardiner (2010), we set c_{s} = Ω_{0} = 10^{3} so that the disk scaleheight is H = Ω_{0}/c_{s} = 1 and one orbital period corresponds to 2π/Ω_{0} in code units. The size of the computational domain, in units of H is and 32 zones per scaleheight are employed. We assume periodic boundary conditions in the azimuthal (y) and vertical (z) directions, while shifted periodicity is imposed at the radial boundary. We carry out two sets of computations, with and without orbital advection, at the resolution of 256 × 256 × 32 zones.
The evolution is characterized by an initial transient phase where all the relevant physical quantities grow rapidly within a few orbits. Subsequently, the flow dynamics becomes nonlinear and the system settles down into a saturated regime accompanied by the formation of typical trailing spiral density waves as shown in Fig. 7 at t = 80π/Ω_{0} (i.e. after 40 revolutions). We note that, owing to the turbulent chaotic behavior of the MRI, a direct comparison between these structures is not really helpful but that one should instead inspect spatially averaged physical quantities. Figure 8 shows the temporal evolution of the volumeintegrated magnetic energy density (top panel) and Maxwell stresses (middle) obtained with and without FARGOMHD for the first 40 orbits. The two computations are in excellent agreement with the absolute value of the (normalized) Maxwell stresses reaching, in both cases, the approximate constant value ≈ 0.01. The time step is plotted in the bottom panel of Fig. 8, where, after the initial transient phase, it is shown that the employment of orbital advection results in a value ≈ 4.3 larger than the standard calculation, in accordance with the results of Stone & Gardiner (2010)
3.4. Diskplanet interaction
Fig. 9 Top: time step as a function of orbits for the diskplanet interaction problem using the standard scheme (black) and FARGO (red). Bottom: change in angular momentum as a function of time. 
Fig. 10 Logarithmic density map for the diskplanet interaction after 100 orbits in the equatorial plane using FARGO (top) and the standard integration scheme (bottom). 
We simulate the interaction of a planet embedded in a viscous global disk as in Kley et al. (2009). The 3D (r,θ,ϕ) computational domain consists of a complete annulus of the protoplanetary disk centered on the star, extending from r_{min} = 0.4 to r_{max} = 2.5 in units of r_{0} = a_{Jup} = 5.2 AU. In the vertical direction, the annulus extends from the disk’s midplane (at θ = 90°) to about 7° (or θ = 83°) above the midplane. The mass of the central star is one solar mass M_{∗} = M_{⊙} and the total disk mass inside [r_{min},r_{max}] is M_{disk} = 0.01 M_{⊙}. For the present study, we use a constant kinematic viscosity coefficient with a value of ν = 10^{15} cm^{2}/s, which corresponds to an αvalue of α = 0.004 at r_{0} for a disk aspect ratio of H = 0.05 rsinθ. The resolution of our simulations is (N_{r},N_{θ},N_{ϕ}) = (256,32,768). At the radial boundaries, the angular velocity is set to the Keplerian values, while we apply reflective, radial boundaryconditions for the remaining variables. In the azimuthal direction we use periodic boundary conditions hold while zerogradient is imposed at the vertical boundaries.
The models are calculated with a locally isothermal configuration where the temperature is constant on cylinders and has the profile T(R) ∝ R^{1} with the cylindrical radius R = rsinθ. This yields a constant ratio of the disk’s vertical height H to the radius R, which was set to 0.05. The initial vertical density stratification is approximately given by a Gaussian (53)Here, the density in the midplane is ρ_{0}(r) ∝ r^{1.5}, which leads to a Σ(r) ∝ r^{ − 1/2} profile of the vertically integrated surface density. The vertical and radial velocities v_{θ} and v_{r} are initialized to zero, while the initial azimuthal velocity v_{ϕ} is given by the equilibrium between gravity, centrifugal acceleration, and the radial pressure gradient. The simulations are performed in a frame of reference corotating with the planet using a secondorder Runge Kutta scheme with linear reconstruction and Courant number C_{a} = 0.25. Following Kley (1998), we adopt a conservative treatment of the Coriolis force in the angular momentum equation.
The planet is described by adding a term to the total gravitational potential acting on the disk as (54)where r_{p} denotes the radius vector of the planet location. We follow Klahr & Kley (2006) in using a cubic planetary potential (55)where r_{sm} = 0.6H. We compute our model using both the standard integration method and the proposed FARGOMHD scheme (without magnetic field). The employment of the FARGO scheme increased the timestep by a factor of about 3.8 (top panel in Fig. 9).
Owing to our choice of boundary conditions and the presence of the planet, we do not fully conserve the total angular momentum, although the angular momentum fluctuations in time, computed with or without orbital advection, are in good agreement (bottom panel in Fig. 9).
Fig. 11 Profiles of surface density (top) and temperature (bottom) as a function of radius for the diskplanet interaction problem. Solid red and black lines refer to computations obtained using FARGO and the standard integration method, respectively. 
In Fig. 10, we show logarithmic maps of density in the equatorial plane after 100 orbits for both of the computations obtained with and without orbital advection. The differences are only marginal. As a more precise test, we also plot the profiles of surface density and temperature as a function of radius, averaged over the azimuthal direction (Fig. 11). The global slope of the profiles as well as the position and width of the gap created by the planet are in excellent agreement between both runs.
3.5. Linear MRI in global disk
We test the FARGOMHD scheme on the linear MRI in global disks, following the setup presented by Flock et al. (2010). This test is very sensitive to the numerical consistency and, at the same time, the MRI growth rates are affected by the numerical dissipation of the scheme.
An analytical description of the linear stage of MRI was given by Balbus & Hawley (1991), while global simulations as well as the nonlinear evolution of the MRI was presented in Hawley & Balbus (1991). The absolute limit of the growth rate for ideal MHD is given for the zero radial wavevectors q_{r} = 0 with the normalized wave vector q_{z}(56)with the Alfvén velocity . The critical mode q_{z} = 0.97 grows exponentially (Ψ = Ψ_{0}e^{γt}) with growth rate γ = 0.75Ω. For global disk models with disk thickness H, the critical wavelength can be rewritten to (57)with the angular frequency Ω = R^{1.5} and the Alfvén velocity V_{A} (see also Eq. (2.3) in Hawley & Balbus 1991).
We use polar cylindrical coordinates (R,φ,Z) with uniform resolution in the region 0 ≤ φ ≤ π/3, − R_{0}/2 ≤ z ≤ R_{0}/2, and R_{0} ≤ R ≤ 4R_{0}, where R_{0} is the unit length. The initial density ρ and pressure p are constant across the entire disk patch with ρ = 1.0, , c_{s} = 0.1V_{φ,0}, and Γ = 1.00001. The gas is initially set up with a Keplerian speed . A uniform vertical magneticfield is placed in the region 2R_{0} ≤ R ≤ 3R_{0}. We choose the strength of the vertical magnetic field to obtain the four fastestgrowing modes, fitting in the domain at R = 2, B_{z} = B_{0}/n where B_{0} = 0.055 and n = 4. We use three different resolutions [R,φ,Z] = [224,168,64] , [112,84,32] , [66,42,16] to measure the capability to resolve the MRI wavelength with 16, 8, and 4 grid cells per vertical scaleheight, respectively. We choose V_{R}(z) = V_{0}sin4z/H as the initial radial velocity for the vertical size H of the box. Boundary conditions are periodic for all variables in the vertical and azimuthal directions, while zero gradient is imposed on all flow quantities at the radial boundaries.
We employ the HLLD Riemann solver Miyoshi & Kusano (2005), piecewise linear reconstruction and the secondorder Runge Kutta timestepping scheme with and without FARGOMHD. The Courant number is set to 0.33 and the spatial reconstruction of the electromotive force at the zone edges is carried out using the approach described in Gardiner & Stone (2005)
Growth rates of the linear MRI mode.
In Table 2, we present the growth of radial magnetic energy over radius during the linear MRI phase using the standard scheme, as well as FARGOMHD, with piecewise linear (Eq. (29)) and piecewise parabolic reconstruction (Eq. (30)). We determine the growth rate from the time derivative of the amplitude maxima for B_{R} in Fourier space at each radius. Computations obtained with the FARGOMHD scheme provide, at the same resolution, a higher growth rate than the standard scheme. This trend has to be attributed to the reduced numerical dissipation as confirmed by the progressive increase in the growth rate when moving from a secondorder to a thirdorder interpolation scheme. The observed speedup is around 4. We also show that we need at least 8 grid cells per wavelength to correctly resolve the growth rate. The thirdorder reconstruction method in combination with FARGOMHD reaches with 16 grid cells per H, growth rates up to the analytical limit of 0.75.
3.6. Turbulent accretion disk
Fig. 12 Threedimensional crosssections of the vertically stratified magnetized accretion disk at t = 100 for the highresolution run showing the surface map of the square of magnetic field B^{2} (left) and density (right). 
As a final application, we consider a turbulent accretion disk in 3D spherical coordinates (r,θ,φ) using the setup BO described by Flock et al. (2011). The initial condition consists of a vertically stratified structure in a Newtonian central potential Φ = −1/r with density (58)where ρ_{0} = 1, R = rsin(θ) is the cylindrical radius, and H/R = c_{0} = 0.07 defines the ratio of scaleheight to radius. The pressure follows locally the isothermal equation of state where the sound speed . The azimuthal velocity is set to (59)The initial velocities V_{r} and V_{θ} are set to a white noise perturbation with amplitude 10^{4}c_{s}. The simulation uses a pure toroidal seed magnetic field with constant plasma β = 2P/B^{2} = 25.
For the sake of comparison, we employ the same computational domain, numerical resolution, and boundary conditions adopted by Flock et al. (2011). We thus have 1 ≤ r ≤ 10 (in AU), π/2 − 0.3 ≤ θ ≤ π/2 + 0.3 (approximately ± 4.3 disk scaleheights), and 0 ≤ φ ≤ 2π. The grid resolution is uniform with N_{r} = 384, N_{θ} = 192, and N_{φ} = 768 zones in the three directions. Buffer zones extend from 1 AU to 2 AU as well as from 9 AU to 10 AU. In the buffer zones, we use a linearly increasing resistivity (up to η = 10^{3}) reaching the boundary. This dampens magnetic field fluctuations and suppresses the interactions with the boundary. Our outflow boundary condition projects the radial gradients in density, pressure, and azimuthal velocity into the radial boundary, and the vertical gradients in density and pressure at the θ boundary. For the first run, we employ the HLLD Riemann solver, piecewise linear reconstruction and 2nd order Runge Kutta time integration. The second run is repeated with the FARGOMHD scheme described in this paper.
As the evolution proceeds, the initial configuration becomes vulnerable to the MRI instability, which quickly leads to a turbulent behavior. Figure 12 shows a slicecut of the square of magnetic field and density after 600 orbits.
A direct comparison between the results obtained with the FARGOMHD scheme and those of Flock et al. (2011) is provided in the three panels of Fig. 13 by plotting relevant volume and timeaveraged quantities. For our analysis, we use the range between 3 AU and 8 AU (in r), which is unaffected by the buffer zones. Temporal averages are taken between 600 and 1200 inner orbits. A detailed description of the measurement can be found in Sect. 2.1 in Flock et al. (2011).
The volumeaveraged values ⟨α_{SS}⟩ of the Shakura and Sunyaev parameter defined by (60)are plotted, as a function of time, in the top panel of Fig. 13. The FARGOMHD run shows a much faster increase in the stresses even within the first 200 inner orbits, at the linear phase of MRI. While the standard computation shows a timeaveraged α_{SS} value of 5 × 10^{3}, the results obtained with FARGOMHD reveal a factor of ≈ 2 increase (α_{SS} = 9.6 × 10^{3}).
The spectra of B_{φ}(m) over azimuthal wavenumber m, plotted in the middle panel of Fig. 13, reveals smaller resolved scales (within a factor of two) of the turbulent magnetic field.
The unstratified global simulations of Sorathia et al. (2012) suggest that the magnetic tilt angle could be a reliable indicator of numerically resolved MRI. We measure the tilt angle for the magnetic field, which is defined by sin2θ_{B} =  B_{r}B_{φ}  /B^{2}, at 4.5 AU and plot the timeaveraged vertical profile of θ_{B} in the bottom panel of Fig. 13. The tilt angle presents the highest values in the coronal region. Employment of the proposed orbitaladvection scheme results overall in a larger tilt angle, of ≈ 10° in the midplane, compared to the standard computation, where ≈ 8° in the midplane. These results suggest that, at the same resolution, the MRI is more accurately resolved when using orbital advection.
With the given grid resolution, chosen to match the results of Flock et al. (2011), we observe a timestep speedup of ~3.75 when employing FARGOMHD. Nonetheless, we note that, according to the guidelines given in Sect. 2.2, this choice may not be an optimal one since the cell aspect ratio is 1:0.67:1.74 (Δr:rΔθ:rΔφ), and larger gains may be obtained by suitably readjusting the number of points in order to have cells with aspect ratios closer to one.
Fig. 13 Top: time history of the accretion stress in the turbulent accretion disk problem. Middle: magnetic field spectra over azimuth. The FARGO method resolves much smaller turbulent scales using the same resolution. Bottom: magnetic tilt angle over height. 
4. Summary
We have presented an orbital advection scheme suitable for the numerical simulations of magnetized, differentially rotating flows in different systems of coordinates. The algorithm has been implemented in the release 4.0 of the PLUTO code for astrophysical fluid dynamics (Mignone et al. 2007, 2012) and shares the same ideas as proposed by Masset (2000), which consist of decomposing, at each timestep, the total velocity into an azimuthally averaged, mean contribution and a residual term. By taking advantage of operator splitting, the two contributions are carried out as a regular update of the standard MHD equations written in the residual velocity and a linear transport step corresponding to a noninteger shift of zones. During the former, any dimensionally unsplit scheme may be employed provided that additional source terms are taken into account and correctly discretized to preserve conservation of both total angular momentum and total energy to machine precision. In both steps, the magnetic field is evolved using the constrained transport formalism to maintain the divergencefree condition.
This approach yields substantially larger time steps whenever the orbital speed exceeds any other characteristic wave signal since the Courant condition depends on the residual velocity rather than the total velocity. Numerical tests confirm that the proposed orbital advection scheme yields results that are equally accurate to and less dissipative than the standard numerical approach at a reduced numerical cost. The overall gain is problemdependent and can be optimized by a suitable choice of grid resolution in the specified geometry: for the selected problems, we observed speedup factors between 4 and 12.
Paralleldomain decomposition can be performed in all three coordinate directions thus allowing the algorithm to be efficiently employed on large numbers of processors on modern parallel computers. The proposed FARGOMHD scheme is particularly suited to largescale global disk simulations that have only recently become amenable to numerical computations on petascale systems.
Acknowledgments
This work has been supported by the PRININAF 2010 grant. We acknowledge the CINECA Award N. HP10CA62TX, 2012 for the availability of highperformance computing resources and support. We gratefully thank the bwGRiD project for the computational resources. A.M. wishes to thank A. Tevzadze for helpful comments on the hydrodynamics vortex problem, G. Mamatsashvili for kindly providing us the reference solution for the MHD shearing wave test, and G. Bodo for valuable comments.
References
 Armitage, P. J. 1998, ApJ, 501, L189 [NASA ADS] [CrossRef] [Google Scholar]
 Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Beckers, J. M. 1992, SIAM J. Num. Anal., 29, 701 [CrossRef] [Google Scholar]
 Beckwith, K., Armitage, P. J., & Simon, J. B. 2011, MNRAS, 416, 361 [NASA ADS] [Google Scholar]
 Bodo, G., Tevzadze, A., Chagelishvili, G., et al. 2007, A&A, 475, 51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Colella, P. 1990, J. Comput. Phys., 87, 171 [NASA ADS] [CrossRef] [Google Scholar]
 Colella, P., & Woodward, P. R. 1984, J. Comput. Phys., 54, 174 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Courant, R., Friedrichs, K., & Lewy, H. 1928, Math. Annal., 100, 32 [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., & Mignone, A. 2010, A&A, 516, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Dzyurkevich, N., Klahr, H., Turner, N. J., & Henning, T. 2011, ApJ, 735, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2005, J. Comput. Phys., 205, 509 [Google Scholar]
 Gardiner, T. A., & Stone, J. M. 2008, J. Comput. Phys., 227, 4123 [NASA ADS] [CrossRef] [Google Scholar]
 Gressel, O., & Ziegler, U. 2007, Comp. Phys. Comm., 176, 652 [Google Scholar]
 Hawley, J. F. 2000, ApJ, 528, 462 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., & Balbus, S. A. 1991, ApJ, 376, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, B. M. 2007, ApJ, 660, 1375 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, B. M., Guan, X., & Gammie, C. F. 2008a, ApJS, 179, 553 [NASA ADS] [CrossRef] [Google Scholar]
 Johnson, B. M., Guan, X., & Gammie, C. F. 2008b, ApJS, 177, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kley, W. 1998, A&A, 338, L37 [Google Scholar]
 Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F. 2000, A&AS, 141, 165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [Google Scholar]
 Mignone, A., Tzeferacos, P., & Bodo, G. 2010, J. Comput. Phys., 229, 5896 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Regev, O., & Umurhan, O. M. 2008, A&A, 481, 21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sorathia, K. A., Reynolds, C. S., Stone, J. M., & Beckwith, K. 2012, ApJ, 749, 189 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137 [NASA ADS] [CrossRef] [Google Scholar]
 Toro, E. 1999, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Applied mechanics: Researchers and students (Springer) [Google Scholar]
 Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1984, SIAM Journal on Scientific and Statistical Computing, 5, 1 [Google Scholar]
Appendix A: Equations in different systems of coordinates
In this section, we give explicit expressions for the MHD Eqs. (8) through (11) written in terms of the residual velocity v′ = v − w for different systems of coordinates. Also, to make notations more compact, we use m = ρv and m′ = ρv′ to denote the total and residual momentum and ℰ′ = −v′ × B will be adopted as a shorthand notation for the (residual) electric field.
A.1. Cartesian coordinates
In cartesian coordinates (x,y,z), we assume, without any loss of generality, that the bulk orbital motion takes place along the y direction, i.e., . The magnitude of w can be defined either analytically or by averaging w in the y direction. In any case, w ≡ w(x,z) so that ∂_{y}w = 0. Writing the equations explicitly (A.1)The source terms on the righthand side of the momentum and energy equations may be written using the reduced forms given in Eqs. (13) and (14). In this case, one has (A.2)However, following the guidelines given in Sect. 2.3.1, a more convenient discretization is given by using Eqs. (21) and (22). This affects only the source terms appearing in both the ycomponent of the momentum equation and in the energy equation which can each now be cast as (A.3)
A.1.1. Shearingbox equations
As a particular case, we briefly review the equations of the shearingbox model introduced in Sect. 3.3. By adopting a noninertial frame that corotates with the disk at orbital frequency Ω_{0}, the momentum and energy Eqs. (3) and (4) become, respectively, (A.4)where is the tidal expansion of the effective gravity while q is the shear parameter (Eq. (47)) and the second term in Eq. (A.4) represents the Coriolis force. The continuity and induction equations retain the same form as the original system in Eqs. (2) and (5).
The derivation of Eq. (A.4) with FARGOMHD is done similarly to the previous section with leading to (A.5)where the source term S_{m′} and S_{E′} can be shown to be, respectively, equal to We note that only the vertical component of gravity is included in the orbital frame (Stone & Gardiner 2010).
A.2. Polar coordinates
In polar cylindrical coordinates (R,φ,z), we consider the bulk orbital velocity to be aligned with the azimuthal direction, that is, where Ω = Ω(R,z) is the angular rotation velocity defined by averaging v_{φ}/R along the azimuthal direction. Writing the equations in components yields (A.8)where, besides the usual divergence operator ∇·( ), we have introduced the “augmented” divergence (A.9)which avoids the curvature source terms in the φcomponent of the residual momentum equation, and is more appropriate for ensuring total angularmomentum conservation.
The source terms in the momentum and energy equations can be written, using the reduced form Eqs. (13) and (14), as (A.10)We note that the total velocity v_{φ} appears in the radial momentum source term. To restore full conservation of total angular momentum and total energy, the formulation given by Eqs. (21) and (22) is more appropriate. This leads to (A.11)
A.3. Spherical coordinates
In spherical coordinates (r,θ,φ) the direction of orbital motion coincides with the azimuthal direction, i.e., where Ω ≡ Ω(r,θ) is total angular rotation velocity computed by averaging v_{φ}/(rsinθ) in the φ direction. Writing Eqs. (8) through (11) in components yields (A.12)where, in analogy with the polar system of coordinates, we have introduced the “augmented” divergence operator (A.13)which is more convenient for expressing conservation of total angular momentum.
The source terms in the momentum and energy equations can be written, using the reduced form Eqs. (13) and (14), as (A.14)
Again, in order to enforce conservation of both total angular momentum and total energy, the formulation given by Eqs. (21) and (22) is the more appropriate to adopt. In this case, one rewrites the φmomentum and energy source terms as (A.15)
All Tables
Average time step for the hydrodynamical vortex test problem at different resolutions using the standard scheme (2nd Col.) and with orbital advection (3rd Col.).
All Figures
Fig. 1 Local potential vorticity for the hydrodynamical vortex test after 20 orbits at the resolution of 1024 × 4096. Top panel: standard integration. Bottom panel: orbital advection. 

In the text 
Fig. 2 Time evolution of the vortensity minimum (top, normalized to its initial value) and relative variation in the total integrated vorticity (bottom) for the hydrodynamical vortex problem. Dotted, dashed, and solid lines refer to computations obtained at the resolutions of 256 × 1024, 512 × 2048, and 1024 × 4096 grid zones, respectively. Black and red colors correspond to the the standard integration algorithm and FARGOMHD, respectively. 

In the text 
Fig. 3 Evolution of the magnetic pressure distribution, ( × 100), for the MHD vortex test computed using the standard MHD update (top row) and with FARGOMHD (bottom row) using 512^{2} zones. From left to right, we show the evolution from time t = 0 to t = 3. 

In the text 
Fig. 4 Top panel: horizontal cuts at y = 0 of the y component of magnetic field at t = 3 for the MHD vortex test in cartesian coordinates using the standard integration method (solid line) and FARGOMHD (plus signs). Bottom panel: timestep variation during the computation. Solid and dashed lines correspond to standard and FARGOMHD computations, respectively. 

In the text 
Fig. 5 Temporal evolution of the volumeaveraged y component of momentum (black) and fractional average energy variation (green) for the MHD vortex with random perturbation using the nonconservative (dashed lines) and conservative (solid line) versions of the orbital advection scheme. 

In the text 
Fig. 6 Azimuthal magneticfield perturbation as a function of time for the compressible MHD shearingwave problem. The exact solution is plotted as a solid black line and the dashed lines correspond to computations obtained at the resolutions of N_{z} = 8 (blue), 16 (green), and 32 (red) zones in the vertical direction. 

In the text 
Fig. 7 Density distribution in the shearingbox computation after 40 rotation periods. Top and bottom panel refers to the model with and without FARGOMHD, respectively. 

In the text 
Fig. 8 Evolution of the MRI in the shearingbox model during the first 40 orbits. From top to bottom, the three panels show, respectively, the volumeintegrated magnetic energydensity, the Maxwell stresses, and the time increment used in the computations without (solid line) and with FARGOMHD (dashed line). 

In the text 
Fig. 9 Top: time step as a function of orbits for the diskplanet interaction problem using the standard scheme (black) and FARGO (red). Bottom: change in angular momentum as a function of time. 

In the text 
Fig. 10 Logarithmic density map for the diskplanet interaction after 100 orbits in the equatorial plane using FARGO (top) and the standard integration scheme (bottom). 

In the text 
Fig. 11 Profiles of surface density (top) and temperature (bottom) as a function of radius for the diskplanet interaction problem. Solid red and black lines refer to computations obtained using FARGO and the standard integration method, respectively. 

In the text 
Fig. 12 Threedimensional crosssections of the vertically stratified magnetized accretion disk at t = 100 for the highresolution run showing the surface map of the square of magnetic field B^{2} (left) and density (right). 

In the text 
Fig. 13 Top: time history of the accretion stress in the turbulent accretion disk problem. Middle: magnetic field spectra over azimuth. The FARGO method resolves much smaller turbulent scales using the same resolution. Bottom: magnetic tilt angle over height. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.