© 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 so-called shearing-box 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 proto-planetary 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 Courant-Friedrichs-Lewy (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 proto-planetary disks where the Keplerian velocity is the dominating flow speed.

In the context of hydrodynamics, Masset (2000) presented the first fast Eulerian-transport 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 co-rotating 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 shearing-box model which is a local approximation describing a small disk patch embodied by a cartesian box co-rotating with the disk at some fiducial radius. Because of the difficulties inherent in simulating an entire disk, local shearing-box models have largely contributed to much of what is presently known about the magneto-rotational 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 magneto-rotational instability, e.g., Armitage (1998), Hawley (2000), Fromang & Nelson (2006) and Regev & Umurhan (2008). Owing to the increased numerical power, high-resolution 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 Godunov-type 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 constrained-transport (CT) scheme, which preserves the divergence-free condition to machine accuracy at all times. The proposed FARGO-MHD 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 time-step limitations imposed by standard explicit time-stepping methods (Sect. 2.1). The new FARGO-MHD 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 angular-momentum 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 FARGO-MHD 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 is the total pressure accounting for thermal (p), and magnetic (B²/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 finite-volume 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 three-dimensional (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 time-stepping scheme. In the PLUTO code, one can take advantage of either the corner-transport upwind (CTU) method of Colella (1990) and Gardiner & Stone (2008) or resort to a Runge-Kutta (RK) discretization, where the spatial gradients are treated as the right-hand 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 6-solve CTU and the second-order Runge-Kutta (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 super-fast, grid-aligned flows, Eq. (7) shows that both CTU and RK2 yield comparable time steps in two-dimensions, while, in three-dimensions, RK2 is expected to give time steps that are approximately twice as large.

Nevertheless, numerical computations of advection-dominated 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 right-hand 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 left-hand 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ⁿ,tⁿ), which corresponds to a non-integer 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 finite-volume or finite-difference 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 two-point 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 point-wise contributions (e.g. gravitational or curvilinear terms) and terms involving the derivatives of the orbital velocity w. Point-wise source terms not involving spatial derivatives are added to the right-hand 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 time-step update. Using, for simplicity, cartesian coordinates and neglecting gravity, a straightforward combination of the equations of continuity and the y-component of momentum yields (20)Although the sum of the last two terms on the right-hand 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 radial-momentum 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 y-momentum upwind fluxes already at disposal during the conservative update. For instance, we update density, (residual) azimuthal-momentum, 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 non-intertial 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ⁿ,t^n + 1]  as q(y,t) = q(y − w(t − tⁿ),tⁿ). 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ⁿ(ξ) = q(ξ,tⁿ) and w does not depend on y. The integrals on the right-hand 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 two-point 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 well-known second-order MUSCL-Hancock 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 δ²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 half-integer 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 divergence-free 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 time-integrated 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 half-increment 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 y-direction, 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 time-integrated 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 time-step gain is problem-dependent 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 time-stepping, one can estimate a time-step 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 sub-grids, which are then solved simultaneously by several processing units. For time-explicit calculations, only boundary data stored in the ghost zones have to be exchanged between neighboring processors. With FARGO-MHD, 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 two-dimensional Keplerian flow

We begin by considering the dynamical evolution of an anti-cyclonic vortex in a two-dimensional (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²), 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 point-mass gravity Φ =  −1/R and fills the computational domain defined by 0.4 ≤ R ≤ 2 and 0 ≤ φ ≤ 2π. With these definitions, the disk scale-height 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 scale-height at R = 1, κ =  −1 is the vortex amplitude, and x = Rcosφ − R₀cosφ₀ and y = Rsinφ − R₀sinφ₀ are the cartesian coordinates measured from the center of the vortex initially located at R₀ = 1 and φ₀ = π/4.

We solve the hydrodynamical equations (no magnetic field) using the second-order Runga-Kutta time-stepping 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.

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 FARGO-MHD, 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, long-lived 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 FARGO-MHD (bottom row) using 5122 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 super-fast sheared flow using cartesian coordinates. The initial condition is similar to the one used by Mignone et al. (2010) and consists of a uniform-density (ρ = 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² = x² + y² 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 y-direction while applying reflecting conditions in the x-directions. 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 field-loop 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² grid zones using the PPM+CTU scheme with Courant number C_a = 0.8.

Results with and without the FARGO-MHD 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 one-dimensional cut of the y-component 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 non-conservative variants of the orbital advection algorithm by plotting the volume-average 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 FARGO-MHD (plus signs). Bottom panel: time-step variation during the computation. Solid and dashed lines correspond to standard and FARGO-MHD computations, respectively.

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

3.3. Shearing-box models

The shearing-box 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 co-rotating with the disk at some fiducial radius R₀ where the orbital frequency is Ω₀. 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 large-scale 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 shearing-box 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 shearing-box 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 shearing-waves in an isothermal medium using the shearing-box 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] ³ 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/Ω₀ 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 Ω₀ = 1.

Fig. 6 Azimuthal magnetic-field perturbation as a function of time for the compressible MHD shearing-wave problem. The exact solution is plotted as a solid black line and the dashed lines correspond to computations obtained at the resolutions of Nz = 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 time-dependent 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 semi-analytical prediction.

3.3.2. Nonlinear MRI

Fig. 7 Density distribution in the shearing-box computation after 40 rotation periods. Top and bottom panel refers to the model with and without FARGO-MHD, respectively.

Fig. 8 Evolution of the MRI in the shearing-box model during the first 40 orbits. From top to bottom, the three panels show, respectively, the volume-integrated magnetic energy-density, the Maxwell stresses, and the time increment used in the computations without (solid line) and with FARGO-MHD (dashed line).

In the next example, we investigate the nonlinear evolution of the magneto-rotational instability (MRI) in the shearing-box model. The initial condition consists of a uniform orbital motion in the y direction of constant density ρ = 1 and zero net-flux 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 = Ω₀ = 10^-3 so that the disk scale-height is H = Ω₀/c_s = 1 and one orbital period corresponds to 2π/Ω₀ in code units. The size of the computational domain, in units of H is and 32 zones per scale-height 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π/Ω₀ (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 volume-integrated magnetic energy density (top panel) and Maxwell stresses (middle) obtained with and without FARGO-MHD 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. Disk-planet interaction

Fig. 9 Top: time step as a function of orbits for the disk-planet 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 disk-planet 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₀ = 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¹⁵ cm²/s, which corresponds to an α-value of α = 0.004 at r₀ 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 boundary-conditions for the remaining variables. In the azimuthal direction we use periodic boundary conditions hold while zero-gradient 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 ρ₀(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 co-rotating with the planet using a second-order 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 FARGO-MHD scheme (without magnetic field). The employment of the FARGO scheme increased the time-step 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 disk-planet 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 FARGO-MHD 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 wave-vectors 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 (Ψ = Ψ₀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₀/2 ≤ z ≤ R₀/2, and R₀ ≤ R ≤ 4R₀, where R₀ 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 magnetic-field is placed in the region 2R₀ ≤ R ≤ 3R₀. We choose the strength of the vertical magnetic field to obtain the four fastest-growing modes, fitting in the domain at R = 2, B_z = B₀/n where B₀ = 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 scale-height, respectively. We choose V_R(z) = V₀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), piece-wise linear reconstruction and the second-order Runge Kutta time-stepping scheme with and without FARGO-MHD. 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)

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 FARGO-MHD, with piece-wise linear (Eq. (29)) and piece-wise 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 FARGO-MHD 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 second-order to a third-order 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 third-order reconstruction method in combination with FARGO-MHD reaches with 16 grid cells per H, growth rates up to the analytical limit of 0.75.

3.6. Turbulent accretion disk

Fig. 12 Three-dimensional cross-sections of the vertically stratified magnetized accretion disk at t = 100 for the high-resolution run showing the surface map of the square of magnetic field B2 (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 ρ₀ = 1, R = rsin(θ) is the cylindrical radius, and H/R = c₀ = 0.07 defines the ratio of scale-height 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^-4c_s. The simulation uses a pure toroidal seed magnetic field with constant plasma β = 2P/B² = 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 scale-heights), 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, piece-wise linear reconstruction and 2nd order Runge Kutta time integration. The second run is repeated with the FARGO-MHD 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 slice-cut of the square of magnetic field and density after 600 orbits.

A direct comparison between the results obtained with the FARGO-MHD scheme and those of Flock et al. (2011) is provided in the three panels of Fig. 13 by plotting relevant volume- and time-averaged 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 volume-averaged 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 FARGO-MHD 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 time-averaged α_SS value of 5 × 10^-3, the results obtained with FARGO-MHD 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_rB_φ | /B², at 4.5 AU and plot the time-averaged 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 orbital-advection 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 time-step speedup of  ~3.75 when employing FARGO-MHD. 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 re-adjusting 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 time-step, 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 non-integer 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 divergence-free 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 problem-dependent 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.

Parallel-domain 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 FARGO-MHD scheme is particularly suited to large-scale global disk simulations that have only recently become amenable to numerical computations on petascale systems.

Acknowledgments

This work has been supported by the PRIN-INAF 2010 grant. We acknowledge the CINECA Award N. HP10CA62TX, 2012 for the availability of high-performance 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

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 short-hand 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 ∂_yw = 0. Writing the equations explicitly (A.1)The source terms on the right-hand 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 y-component of the momentum equation and in the energy equation which can each now be cast as (A.3)

A.1.1. Shearing-box equations

As a particular case, we briefly review the equations of the shearing-box model introduced in Sect. 3.3. By adopting a non-inertial frame that co-rotates with the disk at orbital frequency Ω₀, 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 FARGO-MHD 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 angular-momentum 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

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 FARGO-MHD, 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 FARGO-MHD (bottom row) using 5122 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 FARGO-MHD (plus signs). Bottom panel: time-step variation during the computation. Solid and dashed lines correspond to standard and FARGO-MHD computations, respectively.
In the text

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

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

	Fig. 7 Density distribution in the shearing-box computation after 40 rotation periods. Top and bottom panel refers to the model with and without FARGO-MHD, respectively.
In the text

	Fig. 8 Evolution of the MRI in the shearing-box model during the first 40 orbits. From top to bottom, the three panels show, respectively, the volume-integrated magnetic energy-density, the Maxwell stresses, and the time increment used in the computations without (solid line) and with FARGO-MHD (dashed line).
In the text

	Fig. 9 Top: time step as a function of orbits for the disk-planet 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 disk-planet 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 disk-planet interaction problem. Solid red and black lines refer to computations obtained using FARGO and the standard integration method, respectively.
In the text

	Fig. 12 Three-dimensional cross-sections of the vertically stratified magnetized accretion disk at t = 100 for the high-resolution run showing the surface map of the square of magnetic field B2 (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