RAM: Rapid Advection Algorithm on Arbitrary Meshes

The first step of the RAM method requires the interpolation of the nodes Qⁿ _i to obtain the piecewise functions ${q}_{i}^{n}(x)$. Ideally, this interpolation should be of the order of (or higher than) the one used to solve Equation (6), so errors are dominated by the standard advection step. The RAM method does not rely on any specific choice as long as the interpolation is conservative and does not introduce new extrema. In this work we explore two different conservative piecewise monotonic interpolation methods (i.e., do not introduce new maxima/minima) that are often utilized in hydrodynamics. These are the piecewise linear method (PLM) and piecewise parabolic method (PPM) that we explain below. Throughout this section, all the quantities are evaluated at time t_n , so for the sake of simplicity we omit the superscript n.

Piecewise Linear Method.—This method provides an interpolation that is accurate up to second order (see, e.g., Mignone 2014, Section 3.1). Assume q_i to be of the form

$$\begin{eqnarray}&&{q}_{i}(x)={Q}_{i}+{\sigma }_{i}\left(x-{x}_{i}\right),\end{eqnarray} \tag{ 16 }$$

where σ_i is the slope of the linear function. Monotonicity is preserved by utilizing the slope limiter proposed by van Leer (1974), where

$$\begin{eqnarray}&&{\sigma }_{i}=\left\{\begin{array}{l}0\quad \quad \quad \mathrm{if}\ {\rm{\Delta }}{Q}_{i+\tfrac{1}{2}}{\rm{\Delta }}{Q}_{i-\tfrac{1}{2}}\leqslant 0\\ \\ \displaystyle \frac{2{\rm{\Delta }}{Q}_{i+\tfrac{1}{2}}{\rm{\Delta }}{Q}_{i-\tfrac{1}{2}}}{{\rm{\Delta }}{Q}_{i+\tfrac{1}{2}}+{\rm{\Delta }}{Q}_{i-\tfrac{1}{2}}}\quad \mathrm{otherwise}\ \end{array}\right.,\end{eqnarray} \tag{ 17 }$$

with ${\rm{\Delta }}{Q}_{i+\tfrac{1}{2}}=({Q}_{i+1}-{Q}_{i})/({x}_{i+1}-{x}_{i})$.

Piecewise Parabolic Method.—This method, first introduced by Colella & Woodward (1984), interpolates Q up to third order via quadratic polynomials. The piecewise parabolic function is

$$\begin{eqnarray}\begin{array}{rcl}{q}_{i}(x) & = & {Q}_{L,i}+\left(\displaystyle \frac{x-{x}_{i}}{{\rm{\Delta }}{x}_{i}}\right)\left({Q}_{R,i}-{Q}_{L,i}\right)\\ & & +\,{Q}_{6,i}\left(1-\displaystyle \frac{x-{x}_{i}}{{\rm{\Delta }}{x}_{i}}\right)\left(\displaystyle \frac{x-{x}_{i}}{{\rm{\Delta }}{x}_{i}}\right),\end{array}\end{eqnarray} \tag{ 18 }$$

with ${Q}_{6,i}=6\langle Q{\rangle }_{i}-3\left({Q}_{L,i}+{Q}_{R,i}\right)$. To obtain the quadratic polynomials, we need to calculate 〈Q〉_i (see Section 2) and interpolate 〈Q〉_i to the cell interfaces to obtain $\langle Q{\rangle }_{i+\tfrac{1}{2}}={Q}_{L,i+1}={Q}_{R,i}$ (see Appendix A.1). Furthermore, monotonicity constraints have to be applied (see Appendix A.2). It is worth noticing that the coefficients of the quadratic polynomial depend on 〈Q〉_i , Q_L,i , and Q_R,i ; therefore, the errors of the latter limit the accuracy of the interpolation. In this work we assume 〈Q〉_i ≡ Q_i , which is a second-order approximation of 〈Q〉_i . Therefore, we expect our interpolation errors to be smaller than third order. Higher accuracy could be obtained by improving the calculation of 〈Q〉_i (Colella & Sekora 2008). Note, however, that since PPM conserves 〈Q〉_i locally, a higher-order approximation of 〈Q〉_i may not be conservative in the terms defined in this work, i.e., it does not conserve ∑_i Q_i Δx_i over successive time steps. The conservation will depend on the approximation of the spatial derivative of Q. For example, a fourth-order approximation of 〈Q〉_i = Q_i + 1/24 (Q_i+1 − 2Q_i + Q_i−1) is conservative globally over a periodic domain (i.e., it satisfies Equation (11)). In general, the accuracy of our method can be increased provided that a higher-order approximation of 〈Q〉_i is conserved within a time step (e.g., Felker & Stone 2018).

Finally, we note that for both PLM and PPM monotonicity constraints increase truncation errors locally within nonsmooth regions, which in turn reduces the overall accuracy of the interpolation. An example of reduction of accuracy due to this effect is shown in Section 4.1, for the advection of a square profile.

The volume integral ${M}_{i}^{n+1}$ is obtained from the definition (10) combined with the solution (15),

$$\begin{eqnarray}\begin{array}{rcl}{M}_{i}^{n+1} & = & {\displaystyle \int }_{{x}_{i-\tfrac{1}{2}}}^{{x}_{i+\tfrac{1}{2}}}q(x-{v}_{0}{\rm{\Delta }}t,{t}_{n}){dx}.\end{array}\end{eqnarray} \tag{ 19 }$$

After defining $\tilde{x}\equiv x-{v}_{0}{\rm{\Delta }}t$, we can write Equation (19) as

$$\begin{eqnarray}&&{M}_{i}^{n+1}={\displaystyle \int }_{{\tilde{x}}_{i-\tfrac{1}{2}}}^{{\tilde{x}}_{i+\tfrac{1}{2}}}{q}^{n}(x){dx}.\end{eqnarray} \tag{ 20 }$$

We define positive integers, p(i), such that

$$\begin{eqnarray}&&{x}_{p(i)-\tfrac{1}{2}}\leqslant {\tilde{x}}_{i-\tfrac{1}{2}}\lt {x}_{p(i)+\tfrac{1}{2}},\end{eqnarray} \tag{ 21 }$$

which allows Equation (20) to be evaluated as

$$\begin{eqnarray}&&{M}_{i}^{n+1}=\left\{\begin{array}{l}{\displaystyle \int }_{{\tilde{x}}_{i-\tfrac{1}{2}}}^{{\tilde{x}}_{i+\tfrac{1}{2}}}{q}_{p(i)}^{n}{dx},\quad \mathrm{if}\quad p(i)=p(i+1)\\ {\displaystyle \int }_{{\tilde{x}}_{i-\tfrac{1}{2}}}^{{x}_{p(i)+\tfrac{1}{2}}}{q}_{p(i)}^{n}{dx}+\displaystyle \sum _{m=p(i)+1}^{m\lt p(i+1)}{M}_{m}^{n}\\ +{\displaystyle \int }_{{x}_{p(i+1)-\tfrac{1}{2}}}^{{\tilde{x}}_{i+\tfrac{1}{2}}}{q}_{p(i+1)}^{n}{dx},\quad \mathrm{otherwise}\ \end{array}\right.,\end{eqnarray} \tag{ 22 }$$

where p(N) = p(0) owing to periodicity. Here the sum over M_m ⁿ corresponds to the integral of q_m over an entire cell centered at x_m . Because the reconstruction is conservative, we replace this integral with M_m ⁿ .

Equation (22) is the generalization of standard upwind methods for arbitrary domains of dependence and can be interpreted as a conservative remapping of volume-averaged values (e.g., mass) from coordinates x_i to the new coordinates ${\tilde{x}}_{i}={x}_{i}-{v}_{0}{\rm{\Delta }}t$. Note that since Δt is unbound, ${\tilde{x}}_{i}$ may be nowhere near x_i , and/or a grid cell in x may intersect one, two, or multiple grid cells in $\tilde{x}$ and vice versa. However, this remapping is simplified by the index calculation p(i) (Equation (21) and Section 3.2.1) and the two conditionals in Equation (22). If the time step Δt is bounded (e.g., ${\rm{\Delta }}t\lesssim \min ({\rm{\Delta }}{x}_{i})/{v}_{0}$), Equation (22) reduces to the coordinate remap described in Lufkin & Hawley (1993, see their Equation (58)). In this case, p(i) = i, p(i + 1) = i + 1, and the algorithm in Equation (22) can be cast into the form of a standard upwind advection scheme, similar to the one used to solve Equation (6).

3.2.1. Efficient Calculation of p (i)

If it happens that a coordinate transformation u(x) exists such that ${\rm{\Delta }}u\equiv u({x}_{i+\tfrac{1}{2}})-u({x}_{i-\tfrac{1}{2}})$ is constant, Equation (21) can be written as

$$\begin{eqnarray}&&0\leqslant {\tilde{u}}_{i-\tfrac{1}{2}}-{u}_{p(i)-\tfrac{1}{2}}\lt {\rm{\Delta }}u.\end{eqnarray} \tag{ 23 }$$

In other words, p(i) corresponds to the index of the closest neighbor of the original mesh that is smaller than ${\tilde{u}}_{i-\tfrac{1}{2}}$. Because u is uniformly spaced, this calculation is straightforward,

$$\begin{eqnarray}&&p(i)=\mathrm{int}\left(\displaystyle \frac{{\tilde{u}}_{i}-{u}_{0}}{{\rm{\Delta }}u}\right),\end{eqnarray} \tag{ 24 }$$

where "int(a)" represents the smaller integer that is closer to a. In what follows we introduce a method to build a mesh such that finding this transformation is trivial.

Consider the mesh density function $\psi :[{x}_{\min },{x}_{\max }]\to {{\mathbb{R}}}_{\gt 0}$ such that

$$\begin{eqnarray}&&{\int }_{{x}_{\min }}^{{x}_{\max }}\psi (s)\,{ds}=1.\end{eqnarray} \tag{ 25 }$$

We define $F:[{x}_{\min },{x}_{\min }]\to [0,N-1]$ such that

$$\begin{eqnarray}&&F(x)=(N-1){\int }_{{x}_{\min }}^{x}\psi (s){ds}.\end{eqnarray} \tag{ 26 }$$

The cumulative function F serves as the "index function," as it returns the (continuous) index of a mesh distributed according to the mesh density ψ. This correspondence is evident from the fact that the number of cells per unit space is proportional to dF/dx = ψ. The mesh nodes are defined implicitly by

$$\begin{eqnarray}&&F({x}_{i+\tfrac{1}{2}})-F({x}_{i-\tfrac{1}{2}})=(N-1){\int }_{{x}_{i-\tfrac{1}{2}}}^{{x}_{i+\tfrac{1}{2}}}\psi (s){ds}=1,\end{eqnarray} \tag{ 27 }$$

with ${x}_{-\tfrac{1}{2}}\equiv {x}_{\min }$. Finding analytical expressions for F or its inverse may not be possible, so numerical integration, combined with root-finding methods, has to be used. Finally, the coordinate transformation $u:[{x}_{\min },{x}_{\max }]\to [0,1]$ such that

$$\begin{eqnarray}&&u(x)={\int }_{{x}_{\min }}^{x}\psi (s){ds}\end{eqnarray} \tag{ 28 }$$

ensures that, by definition, Δu is constant, and Equation (24) evaluates to

$$\begin{eqnarray}&&p(i)=\mathrm{int}[F({\tilde{u}}_{i})],\end{eqnarray} \tag{ 29 }$$

which allows us to speed up index calculations needed in step 2 of RAM significantly. It is worth noticing that Hugger (1993) generalizes a similar method to multiple dimensions.

Once ${M}_{i}^{n+1}$ is obtained, the third step consists of using Equation (10) to obtain the advected profile ${Q}_{i}^{n+1}$ as

$$\begin{eqnarray}&&{Q}_{i}^{n+1}=\displaystyle \frac{{M}_{i}^{n+1}}{{\rm{\Delta }}{x}_{i}}.\end{eqnarray} \tag{ 30 }$$

In Figure 1 we advect to the right a sinusoidal profile, Q(x, t_n ), sampled on a discrete nonuniform grid over a length v₀Δt. To advect it, we follow the three steps described previously. The first step (top panel) consists in calculating q_i ⁿ , which, for simplicity, is done using PLM. These linear functions are plotted with red lines. To explain the process in more detail, we now focus on the interval $[{x}_{4-\tfrac{1}{2}},{x}_{4+\tfrac{1}{2}}]$ only, but a similar procedure has to be done for all the intervals. The goal is to find ${Q}_{4}^{n+1}$ at x₄ and time t_n+1 that results from advecting qⁿ . For that we follow step 2 (middle panel). We transform x₄ to ${\tilde{x}}_{4}$ and find all p(i) that are within $[{\tilde{x}}_{4-\tfrac{1}{2}},{\tilde{x}}_{4+\tfrac{1}{2}}]$. This step can be done efficiently if the mesh follows a mesh density function ψ as explained in Section 3.2.1. For this particular example, the specific indices needed to evaluate Equation (22) in this interval are p(4), p(5). To satisfy the inequality (21), p(4) = 1 and p(5) = 3, and Equation (22) becomes

$$\begin{eqnarray}&&{M}_{4}^{n+1}={\int }_{{\tilde{x}}_{4-\tfrac{1}{2}}}^{{x}_{1+\tfrac{1}{2}}}{q}_{1}^{n}\,{dx}+{M}_{2}^{n}+{\int }_{{x}_{3-\tfrac{1}{2}}}^{{\tilde{x}}_{4+\tfrac{1}{2}}}{q}_{3}^{n}\,{dx}.\end{eqnarray} \tag{ 31 }$$

Thus, to calculate ${M}_{4}^{n+1}$, we have to integrate q₁ ⁿ , q₂ ⁿ , and q₃ ⁿ in the interval $[{\tilde{x}}_{4-\tfrac{1}{2}},{\tilde{x}}_{4+\tfrac{1}{2}}]$. Here the integral of q₂ ⁿ is replaced by M₂ ⁿ simply because the step of building q₂ is conservative. Finally, in step 3 (bottom panel) we take ${M}_{4}^{n+1}$ and use Equation (30) to calculate ${Q}_{4}^{n+1}$ (depicted with an open triangle) at x₄ as ${Q}_{4}^{n+1}={M}_{4}^{n+1}/{\rm{\Delta }}{x}_{4}$.

RAM can be implemented in two- and three-dimensional hydrodynamics solvers through dimensional splitting (Stone & Norman 1992; Benítez-Llambay & Masset 2016, Section 3.4) and used to advect along the direction for which the bulk velocity is larger. In this case, an additional limitation to the CFL condition due to relative shear needs to be implemented (see Benítez-Llambay & Masset 2016, Section 3.5.2). In Section 4.2, we show an example of a two-dimensional polar disk in which RAM was implemented through dimensional splitting, using it to solve advection along the azimuthal direction, for which the Keplerian speed dominates the dynamics. RAM is also agnostic to the discretization scheme. It can be used for orbital advection in both finite-difference and finite-volume methods, as done previously with the uniform-mesh FARGO method (Johnson et al. 2008; Stone & Gardiner 2010; Mignone et al. 2012; Lyra et al. 2017). We finally note that if higher-order time integration is needed on multidimensional problems, an extra source term needs to be considered in the hydrodynamics equations (see Shariff & Wray 2018, Equation (33)).

To compare the errors and convergence rates for the different advection methods discussed above, we study the one-dimensional advection of smooth and sharp profiles on uniform and nonuniform meshes. In this test, we consider the transport step of FARGO3D only (see Benítez-Llambay & Masset 2016, for details). This test highlights the errors associated with advection that a general hydrodynamics integration has for different advection algorithms. The periodic domain spans the range [−π, π]. For the smooth profile, the initial condition is

$$\begin{eqnarray}&&Q(x)=\displaystyle \frac{2}{\pi }{e}^{-4{x}^{2}/\pi },\end{eqnarray} \tag{ 32 }$$

whereas for the square profile it is

$$\begin{eqnarray}&&Q(x)=\displaystyle \frac{1}{\pi }\left\{\begin{array}{l}0.75\quad | x| \leqslant \pi /2\\ 0.25\quad \mathrm{otherwise}\end{array}\right..\end{eqnarray} \tag{ 33 }$$

While Equation (32) is not strictly periodic within the domain, the discontinuity produced at the edges in the initial condition is very small and does not introduce any significant error. In the two cases, the advection velocity is v = π, and it is used as the bulk velocity v₀ for both FARGO and RAM. We integrate the system for a time t = 2, so the profile moves a full period, returning to its initial position. To test convergence with resolution, we double the number of grid points from 32 up to 512 and estimate truncation errors, e₁, through the L1-norm,

$$\begin{eqnarray}&&{e}_{1}=\displaystyle \frac{{\sum }_{i=0}^{N-1}| {Q}_{i}-{Q}_{i}^{\mathrm{ref}}| {\rm{\Delta }}{x}_{i}}{{\sum }_{i=0}^{N-1}{\rm{\Delta }}{x}_{i}},\end{eqnarray} \tag{ 34 }$$

where Q^ref is a reference solution, which we take to be equal to the initial condition. This is because after one period the exact solution corresponds to the same initial profile (see Equation (7)). We consider two cases, one in which the mesh is uniformly distributed and another where the mesh is distributed according to the periodic bump mesh density function described in Appendix B, with parameters a = π/2, b − a = 0.1, c = 5 (see also Section 3.2.1). This density function enables a smooth transition between two refinement levels, each with uniform resolution. The time step for the integration is the standard CFL condition

$$\begin{eqnarray}&&{\rm{\Delta }}t=\beta \,\min \left(\displaystyle \frac{{\rm{\Delta }}{x}_{i}}{{v}_{0}}\right),\end{eqnarray} \tag{ 35 }$$

with β the CFL factor, which is usually smaller than 1 because of stability constraints. We arbitrarily choose β = 0.5 for STD ⁷ and β = 5.5 for FARGO and RAM, which is of the order of the upper bound for the typical speed-up obtained with FARGO in protoplanetary disk simulations.

In Figure 2 we summarize the results obtained for this test. We find that when advecting a smooth profile truncation errors for all the methods are second order, as they decrease quadratically with N (e₁ ∝ N⁻²), whereas for the sharp profile the errors scale ∝N^−3/4, consistent with a first-order algorithm. For both uniform and nonuniform meshes, the absolute errors of RAM with PLM are slightly larger than errors with PPM, as expected. For the uniform mesh, errors using FARGO or RAM with PPM are identical and decrease with increasing resolution slightly faster than the errors using STD or RAM with PLM; however, the difference in slope is very small. The errors with RAM (PPM) are the truncation error of the PPM interpolation. On the other hand, the errors using FARGO correspond to an additional advection with the residual velocity ${v}_{{\rm{shift}}}^{R}={v}_{0}-N{\rm{\Delta }}x/{\rm{\Delta }}t$ (because the time step does not allow the profile to be shifted an integer number of cells only). This residual step is solved with the STD advection scheme with PPM reconstruction (see Equation (69) in Benítez-Llambay & Masset 2016). Thus, it is not surprising that FARGO and RAM with PPM produce the same errors.

Zoom In Zoom Out Reset image size

Finally, the large difference between the absolute errors of STD and RAM/FARGO is due to the significant difference between the time steps chosen for each method. As discussed in Section 1, truncation errors scale with the time step and the bulk velocity v₀. This highlights the importance of having a method that allows advection using larger time steps on both uniform and nonuniform domains.

This first test demonstrates the validity of our formulation of RAM. As expected, the precision and convergence properties of RAM are controlled by the interpolation step. It also shows that RAM is equivalent to FARGO on uniform meshes, while being its generalization on nonuniform grids.

In this section, we demonstrate that RAM, combined with the mesh density function, can be used to accelerate hydrodynamics calculations of differentially rotating fluids and decrease numerical diffusion on nonuniform grids. With this purpose in mind, we set up a two-dimensional viscous disk orbiting a central star of mass M_⋆ = M₀. The disk model ⁸ corresponds to a steady-state α-disk (Shakura & Sunyaev 1973), in which the gas surface density, Σ, and temperature, T, are power laws in radius with exponents –1/2 and –1, respectively. The disk aspect ratio and viscosity parameters are set to h = 0.05 and α = 10⁻³. We use polar coordinates, where the disk domain extends radially and azimuthally as $r/{r}_{0}\in \left[0.4,2.1\right]$, $\varphi \in \left[-\pi ,\pi \right]$, respectively. The disk is assumed to be locally isothermal, with a sound speed c_s = hv_K, where ${v}_{{\rm{K}}}=\sqrt{{{GM}}_{\star }/r}$ is the Keplerian speed and G is the gravitational constant. We include an embedded planet of mass m _p = 10⁻³ M_⋆ orbiting the central star on a circular orbit of radius r _p /r₀ = 1 and modeled as a Plummer potential with small softening length 3 × 10⁻³ r₀. The planet is inserted at t = 0, and the calculations are stopped after one orbital period of the planet, at $t=2\pi {{\rm{\Omega }}}_{0}^{-1}$, with ${{\rm{\Omega }}}_{0}=\sqrt{{{GM}}_{0}/{r}_{0}^{3}}$ the Keplerian frequency at r₀. The hydrodynamics equations are solved in a frame that corotates with the planet, where the planet is fixed at an azimuth φ _p = 0. We choose a unit system such that G = r₀ = M₀ = 1. The intent of this test is not to simulate planet−disk interaction realistically but to demonstrate the efficiency of RAM for an astrophysical scenario with high Mach number bulk velocities, high amplitude, and small-scale flow features. We compare the computational cost of RAM with traditional advection schemes.

We run simulations with 45, 90, and 180 cells within the Hill radius ${R}_{{\rm{H}}}={r}_{p}\sqrt[3]{{{Gm}}_{p}/(3{M}_{\star })}$ on both azimuthally uniform and nonuniform meshes. To conserve resolution per scale height across the disk domain, meshes that are evenly spaced in azimuth are also radially uniform in logarithmic space, which is a common choice for the study of planet−disk interaction. Runs with nonuniform spacing are distributed azimuthally according to the periodic bump mesh density function described in Appendix B, with parameters a = 0.1387, b = 0.2773, and c = 15. In the radial direction, we use the modified bump of Appendix B.2, centered at s₀ = 1, with the same a and b parameters, but with c = 2.5. In Figure 3, we depict the structured nonuniform mesh generated with these mesh density functions. The white region corresponds to the resolution transition controlled by the difference b − a. This mesh increases resolution close to the planet in the radial and azimuthal directions. In addition, because of the logarithmic nature of the radial mesh, the cell size decreases inward. When the azimuthal spacing is uniform, we utilize both FARGO and STD as advection methods (we already showed in Section 4.1 that RAM is equivalent to FARGO on uniform domains). Otherwise, when the azimuthal spacing is nonuniform, advection is solved with STD and RAM. In all cases, we use the faster PLM reconstruction method. As demonstrated in Section 4.1, results with PPM are similar but with lower truncation errors. Due to the small azimuthal size of the cells in the innermost region of the disk and the large bulk speed, a RAM-like method is expected to produce a large speed-up.

Zoom In Zoom Out Reset image size

4.2.1. Theoretical Speed-up of RAM

We analyze the main terms in the CFL stability condition to obtain a first estimate of the expected computational gain using RAM for this particular problem. The time step is computed following Section 3.5 in Benítez-Llambay & Masset (2016). We also account for the strong gravitational acceleration produced by the planet by assuming that rotational balance can be reached quickly after the planet is inserted in the disk. In this way, Δt may be limited by the additional local Keplerian speed induced within the Hill radius. This estimation corresponds to a lower bound for Δt because, in reality, rotational equilibrium may take longer to be established fully or pressure support can make the disk rotate slower. The expression we use to calculate Δt is

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}t=\beta \min \left[\left({c}_{{\rm{s}}}+{v}_{{\rm{p}}}\right)\max (\xi ,\eta )+\xi | {v}_{{\rm{K}}}-{v}_{\mathrm{frame}}| \right.\\ {\left.+4\nu \max ({\xi }^{2},{\eta }^{2})\right]}^{-1},\end{array}\end{eqnarray} \tag{ 36 }$$

where ξ = 1/(rΔφ), η = 1/Δr, and v_frame = Ω_K(r_p)r is the azimuthal velocity of the rotating frame. Following the FARGO3D implementation, we define the CFL parameter β = 0.44. The rotational velocity induced by the planet, v_p, is given by

$$\begin{eqnarray}&&{v}_{p}=\sqrt{\frac{{{Gm}}_{{\rm{p}}}}{s}}{e}^{-\tfrac{{s}^{2}}{2{\sigma }^{2}}},\end{eqnarray} \tag{ 37 }$$

with s = ∣r − r_p∣ and σ = R_H/2. Here the Keplerian speed has been multiplied by an exponential cutoff to make the effect of rotation local.

Equation (36) includes the contribution of the bulk velocity of the disk, v_K − v_frame, and the velocity induced by the planet, v_p. As such, for a given mesh, it estimates Δt for STD. To obtain the time step allowed by RAM or FARGO, we set v_K − v_frame = 0. In addition, the maximum speed-up possible by RAM can be calculated after setting v_p = 0. In Figure 4, we show the results of Equation (36) when applied to nonuniform (left panels) and uniform (right panels) meshes with 180 cells/R_H. For each mesh, the top and bottom panels show the allowed time step in a region that includes both the planet and the inner boundary for STD and RAM or FARGO, respectively. Since c_s increases inward while the cell size decreases (see Figure 3), Δt decreases toward the inner boundary. In addition, because in the nonuniform mesh there is an azimuthal band in which the azimuthal size also decreases with decreasing r, Δt can become even more restrictive. This effect is particularly strong for STD, which also has the contribution of the bulk velocity, v_K − v_frame, that increases rapidly for decreasing r. Another remarkable feature is the small Δt produced owing to the high rotational velocity induced by the planet within a high-resolution region around (x, y) = (r₀, 0). Whether or not this velocity dominates the time step of the simulation depends exclusively on the location of the inner boundary (the closer it is to the central star, the smaller the time step). When using RAM, the absence of the term proportional to the bulk velocity allows Δt to be larger almost everywhere, except at the planet's orbit, where the bulk velocity is, by construction, zero. Therefore, in the case of RAM, the terms that limit Δt are those associated with c_s and v_p. The former may limit Δt owing to the double effect of c_s being large inward and the cell size being also small there, while the latter may limit Δt owing to the small cell size and large velocity induced close to the planet. In the case in which the mesh is uniform, we also calculate Δt for FARGO and STD. The time step allowed globally by this mesh is almost indistinguishable from the one obtained for the nonuniform case. The same dependencies of Δt with r are observed. It is also worth noticing that, for both meshes, viscosity could limit Δt strongly if either α or the resolution is further increased.

Zoom In Zoom Out Reset image size

We now report the maximum allowed Δt obtained for RAM and STD with/without a planet on the nonuniform mesh. For STD with a planet, the maximum time step is ${\rm{\Delta }}{t}_{\mathrm{STD}}\,=5.7\times {10}^{-5}{{\rm{\Omega }}}_{0}^{-1}$, while for RAM we obtain ${\rm{\Delta }}{t}_{\mathrm{RAM}}\,=2.9\times {10}^{-4}{{\rm{\Omega }}}_{0}^{-1}$. This shows that RAM allows a time step that is ∼5 times larger than STD. If we do not consider the planet in this calculation, for STD we obtain ${\rm{\Delta }}{t}_{\mathrm{STD},\mathrm{np}}=5.7\times {10}^{-5}{{\rm{\Omega }}}_{0}^{-1}$ and ${\rm{\Delta }}{t}_{\mathrm{RAM},\mathrm{np}}=7.6\times {10}^{-4}{{\rm{\Omega }}}_{0}^{-1}$ for RAM, which is ∼13 times larger than that of STD.

Finally, when considering a uniform mesh, producing the same resolution of 180 cells/R_H requires 16 times more cells than in the nonuniform case, using eight and two times more cells in azimuth and radius, respectively. Given that RAM and FARGO produce the same time step, RAM on a nonuniform mesh is 16 times faster than a FARGO run on a uniform mesh. When comparing RAM on nonuniform meshes with respect to STD on a uniform one, RAM can be ∼80 − 208 times faster than STD on a corotating frame, and even faster in noncorotating frames.

4.2.2. Simulation Results

On large scales, the features obtained with different methods/resolutions are indistinguishable from one another. We highlight these features in Figure 5, which shows the surface density, Σ, for the run advected with RAM and 180 cells/R_H. In Figure 6, we show snapshots Σ for the runs with 45 cells/R_H (left column) and 180 cells/R_H (middle column), for RAM (top panels), FARGO (middle panels), and STD in a nonuniform domain (bottom panels) after one orbit of the planet. In each panel, we zoom in on the circumplanetary region where spirals develop. In addition, in the right panels we plot cuts along the black line drawn in the zoom-in for different resolutions (low/high resolution with blue/orange curves).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is when we look within the circumplanetary region that significant differences are observed for runs with low and high resolution. The ambient density within the Hill radius and the density contrast of the spirals are most impacted. Low-resolution runs are affected by larger truncation errors and numerical diffusion. The effect of these errors is to increase the ambient density within R_H and smooth out the spirals. This can be recognized when comparing low-resolution runs with FARGO, RAM, and STD with both low and high resolution (blue curves in the right panels). We note that the high numerical diffusion of STD within the circumplanetary region is not due to the bulk velocity because this velocity is small in the corotating frame used. Hence, we conclude that the smaller time steps of STD and the corresponding larger number of interpolations/reconstructions required to advance the solution up to the final integration time are the main source of numerical diffusion close to the planet. It is worth noting that RAM on a nonuniform mesh reproduces the sharp features obtained with FARGO. However, the nonuniform mesh utilizes 16 times fewer grid points than FARGO, significantly reducing the computational cost while maintaining the same level of accuracy.

We now report measurements of the time step required for each method in different meshes for the high-resolution run with 180 cells/R_H only, but similar results were obtained for the other runs with lower resolution. We find that ${\rm{\Delta }}{t}_{\mathrm{STD}}=5.7\times {10}^{-5}{{\rm{\Omega }}}_{0}^{-1}$ and ${\rm{\Delta }}{t}_{\mathrm{RAM}}=3.87\times {10}^{-4}{{\rm{\Omega }}}_{0}^{-1}$. Despite the simplified modeling, these values are in excellent agreement with the estimations described in Section 4.2.1. Ideally, the relative speed-up between methods is given by the ratio Δt_RAM/Δt_STD ≃ 6.9; however, RAM introduces an overhead with respect to STD owing to the velocity splitting (see Section 1) and the additional calculations required to shift the profile (see Section 3). Using fast index calculation (see Section 3.2.1), in our tests we found that this overhead is ∼20%. Therefore, the effective speed-up is Δt_RAM/Δt_STD × 0.8 ≃ 5.5. Note that a similar overhead is introduced with FARGO, but it requires eight times more cells in azimuth to match the resolution of 180 cells/R_H. Since we find Δt_RAM ≃ Δt_FARGO, we conclude that FARGO is eight times more expensive than RAM on nonuniform meshes provided that the radial spacing used for both methods is identical. However, in uniform meshes like the one presented in Section 4.2, where the number of cells is two and eight times more in radius and azimuth, respectively, FARGO is 16 times more expensive than RAM on the nonuniform mesh.

The problem studied in this section illustrates the main advantages of RAM to study the dynamics close to a planet at a lower computational cost. Two main results can be highlighted from our analysis. First, the standard approach using FARGO to simulate the problem on an azimuthally uniform mesh is expensive. Interestingly, a cheaper method consists in utilizing STD on a nonuniform mesh (in both radius and azimuth). Second, our results clearly show the advantage of using nonuniform meshes combined with RAM, which has smaller truncation errors and is significantly cheaper than STD.

Following Colella & Woodward (1984), $\langle Q{\rangle }_{i+\tfrac{1}{2}}$ are calculated as

$$\begin{eqnarray}\begin{array}{rcl}\langle Q{\rangle }_{i+\tfrac{1}{2}} & = & \langle Q{\rangle }_{i}+\displaystyle \frac{{\rm{\Delta }}{x}_{i}}{{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}\left(\langle Q{\rangle }_{i+1}-\langle Q{\rangle }_{i}\right)\\ & & +\displaystyle \frac{\left({\tilde{\omega }}_{1}(\langle Q{\rangle }_{i+1}-\langle Q{\rangle }_{i})+{\tilde{\omega }}_{2}\delta \langle Q{\rangle }_{i+1}+{\tilde{\omega }}_{3}\delta \langle Q{\rangle }_{i}\right)}{\displaystyle \sum _{k=-1}^{2}{\rm{\Delta }}{x}_{i+k}},\end{array}\end{eqnarray} \tag{ A1 }$$

with

$$\begin{eqnarray}&&\begin{array}{c}{\tilde{\omega }}_{1}=\displaystyle \frac{2{\rm{\Delta }}{x}_{i+1}{\rm{\Delta }}{x}_{i}}{{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}\left[\displaystyle \frac{{\rm{\Delta }}{x}_{i-1}+{\rm{\Delta }}{x}_{i}}{2{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}-\displaystyle \frac{{\rm{\Delta }}{x}_{i+2}{\rm{\Delta }}{x}_{i+1}}{2{\rm{\Delta }}{x}_{i+1}+{\rm{\Delta }}{x}_{i}}\right]\\ {\tilde{\omega }}_{2}=-{\rm{\Delta }}{x}_{i}\displaystyle \frac{{\rm{\Delta }}{x}_{i-1}+{\rm{\Delta }}{x}_{i}}{2{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}\\ {\tilde{\omega }}_{3}={\rm{\Delta }}{x}_{i+1}\displaystyle \frac{{\rm{\Delta }}{x}_{i+1}+{\rm{\Delta }}{x}_{i+2}}{{\rm{\Delta }}{x}_{i}+2{\rm{\Delta }}{x}_{i+1}}\end{array}\end{eqnarray} \tag{ A2 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\delta \langle Q{\rangle }_{i} & = & \displaystyle \frac{{\rm{\Delta }}{x}_{i}}{{\rm{\Delta }}{x}_{i-1}+{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}\\ & & \times \,\left(\displaystyle \frac{2{\rm{\Delta }}{x}_{i-1}+{\rm{\Delta }}{x}_{i}}{{\rm{\Delta }}{x}_{i}+{\rm{\Delta }}{x}_{i+1}}\left(\langle Q{\rangle }_{i+1}-\langle Q{\rangle }_{i}\right)\right.\\ & & +\,\left.\displaystyle \frac{{\rm{\Delta }}{x}_{i}+2{\rm{\Delta }}{x}_{i+1}}{{\rm{\Delta }}{x}_{i-1}+{\rm{\Delta }}{x}_{i}}\left(\langle Q{\rangle }_{i}-\langle Q{\rangle }_{i-1}\right)\right).\end{array}\end{eqnarray} \tag{ A3 }$$

To apply the monotonicity constraint, first δ〈Q〉_i is modified in the calculation $\langle Q{\rangle }_{i+\tfrac{1}{2}}$ as

$$\begin{eqnarray*}&&\begin{array}{c}\delta \langle Q{\rangle }_{i}=\left\{\begin{array}{c}{\rm{s}}{\rm{g}}{\rm{n}}(\delta \langle Q\rangle i)min\left(|\delta \langle Q{\rangle }_{i}|,2|\langle Q{\rangle }_{i+1}-\langle Q{\rangle }_{i}|,2|\langle Q{\rangle }_{i}-\langle Q{\rangle }_{i-1}|\right)\,{\rm{i}}{\rm{f}}\left(\langle Q{\rangle }_{i+1}-\langle Q{\rangle }_{i}\right)\left(\langle Q{\rangle }_{i}-\langle Q{\rangle }_{i-1}\right)\gt 0\\ 0\,{\rm{o}}{\rm{t}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{w}}{\rm{i}}{\rm{s}}{\rm{e}}\end{array}\right..\end{array}\end{eqnarray*}$$

With updated values of $\langle Q{\rangle }_{i+\tfrac{1}{2}}$, the coefficients Q_R,i and Q_L,i are obtained as

$$\begin{eqnarray}&&{Q}_{R,i}={Q}_{L,i}=\langle Q{\rangle }_{i}\quad \quad \mathrm{if}\,\left({Q}_{R,i}-\langle Q{\rangle }_{i}\right)\left(\langle Q{\rangle }_{i}-{Q}_{L,i}\right)\leqslant 0\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{Q}_{L,i}=3\langle Q{\rangle }_{i}-2{Q}_{R,i}\,\quad \mathrm{if}\,\left({Q}_{R,i}-{Q}_{L,i}\right){Q}_{a,i}\gt \frac{1}{6}{\left({Q}_{R,i}-{Q}_{L,i}\right)}^{2}\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&{Q}_{R,i}=3\langle Q{\rangle }_{i}-2{Q}_{L,i}\,\,\,{\rm{if}}-\frac{1}{6}{\left({Q}_{R,i}-{Q}_{L,i}\right)}^{2}\gt \left({Q}_{R,i}-{Q}_{L,i}\right){Q}_{a,i},\end{eqnarray} \tag{ A6 }$$

with ${Q}_{a,i}=\langle Q{\rangle }_{i}-\tfrac{1}{2}\left({Q}_{L,i}+{Q}_{R,i}\right)$.

In this work, we restrict our analysis to functions F with known analytical expressions for the integrals and use the bisection method to solve Equation (27) or the expression of the inverse when it is easy to find.

We define the function ψ(s) as

$$\begin{eqnarray}\psi (s)=\left\{\begin{array}{lr}1+c & | s| \leqslant a\\ 1+c{\cos }^{2}\left(\displaystyle \frac{\pi }{2}\displaystyle \frac{| s| -a}{b-a}\right) & a\lt | s| \lt b\\ 1 & | s| \geqslant b\end{array}\right.,\end{eqnarray} \tag{ B1 }$$

with ${x}_{\min }\leqslant s\leqslant {x}_{\max }$. To satisfy the condition given by Equation (25), we normalize ψ(s) by I_ψ , which is obtained after integrating the density function over the domain $\left[{x}_{\max },{x}_{\min }\right]$, that is,

$$\begin{eqnarray}&&{I}_{\psi }={\displaystyle \int }_{{x}_{\min }}^{{x}_{\max }}\psi (s){ds}=({x}_{\max }-{x}_{\min })+c\left(a+b\right).\end{eqnarray} \tag{ B2 }$$

Thus, the desired mesh density function corresponds to

$$\begin{eqnarray}&&\bar{\psi }={I}_{\psi }^{-1}\psi .\end{eqnarray} \tag{ B3 }$$

This function generates a mesh with two levels of uniform resolution and includes a smooth transition between the two levels. The value of the parameter a sets the extent of the highest-resolution level. The parameter b sets the location of the transition between two levels, and therefore ∣b∣ > ∣a∣. The extent of the transition is given by b − a. The value of the parameter c sets the relative resolution between the two levels.

For this case, the cumulative $F(x)={\int }_{-\pi }^{x}\bar{\psi }(s){ds}$, which corresponds to the definition of the coordinate u(x), is

$$\begin{eqnarray}\begin{array}{rcl}F(x) & = & \left\{\begin{array}{ll}f(x)-f(-\pi ) & \quad x\leqslant -b\\ {g}^{+}(x)-{g}^{+}(-b)+F(-b) & \quad -b\lt x\leqslant -a\\ h(x)-h(-a)+F(-a) & \quad -a\lt x\leqslant a\\ {g}^{-}(x)-{g}^{-}(a)+F(a) & \quad a\lt x\leqslant b\\ f(x)-f(b)+F(b) & \quad x\lt b\end{array}\right.,\end{array}\end{eqnarray} \tag{ B4 }$$

with the functions f, g^±, h defined as

$$\begin{eqnarray}&&f(x)=\displaystyle \frac{x}{I},\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&{g}^{\pm }(x)=\displaystyle \frac{\left(2+c\right)}{2}f(x)-\displaystyle \frac{c\left(b-a\right)}{2\pi I}\sin \left(\pi \displaystyle \frac{x\pm a}{a-b}\right),\end{eqnarray} \tag{ B6 }$$

$$\begin{eqnarray}&&h(x)=\left(1+c\right)f(x),\end{eqnarray} \tag{ B7 }$$

In Figure 7 we plot an example of a mesh density function with $s\in \left[-\pi ,\pi \right]$. The parameters are defined to provide a ∼5 × resolution contrast between the two levels. In the right panel of Figure 7 we include the cumulative distribution for a mesh with N = 64 points. Finally, we note that the piecewise density function proposed here can be modified to allow for arbitrary levels of resolution with continuous transitions.

Zoom In Zoom Out Reset image size

In Section 4.2 we utilize a mesh density function in the radial direction centered at s₀ = 1 and defined as

$$\begin{eqnarray}\psi (s)=\left\{\begin{array}{lr}1/s+c & | s-{s}_{0}| \leqslant a\\ 1/s+c{\cos }^{2}\left(\displaystyle \frac{\pi }{2}\displaystyle \frac{| s-{s}_{0}| -a}{b-a}\right) & a\lt | s-{s}_{0}| \lt b\\ 1/s & | s-{s}_{0}| \geqslant b\end{array}\right.,\end{eqnarray} \tag{ B8 }$$

with ${I}_{\psi }=\mathrm{log}({x}_{\max })-\mathrm{log}({x}_{\min })+c(a+b)$ (see Equation (B2)). Note that when the parameter c = 0 we obtain a log-uniform mesh density.

Consider a 1D pressure gradient along the y-direction. The velocity is updated as

$$\begin{eqnarray}&&{v}_{{ij}-\tfrac{1}{2}k}^{n+s}={v}_{{ij}-\tfrac{1}{2}k}^{n}-\displaystyle \frac{{\rm{\Delta }}t}{{\left\langle {\rho }^{n}\right\rangle }_{{ij}-\tfrac{1}{2}k}}\displaystyle \frac{{P}_{{ijk}}^{n}-{P}_{{ij}-1k}^{n}}{{y}_{j}-{y}_{j-1}},\end{eqnarray} \tag{ C3 }$$

where n+s indicates a time update from step n and ${\left\langle {\rho }^{n}\right\rangle }_{{ij}-\tfrac{1}{2}k}$ is a spatial averaged density (see Benítez-Llambay & Masset 2016). It can be shown that momentum is conserved in uniform and nonuniform grids if the average density is defined as

$$\begin{eqnarray}&&{\left\langle {\rho }^{n}\right\rangle }_{{ij}-\tfrac{1}{2}k}=\displaystyle \frac{1}{2}\displaystyle \frac{{\rho }_{{ijk}}^{n}\left({y}_{j+\tfrac{1}{2}}-{y}_{j-\tfrac{1}{2}}\right)+{\rho }_{{ij}-1k}^{n}\left({y}_{j-\tfrac{1}{2}}-{y}_{j-\tfrac{3}{2}}\right)}{{y}_{j}-{y}_{j-1}}.\end{eqnarray} \tag{ C4 }$$

The same holds for cylindrical and spherical coordinates. Note that the same update needs to be done in substep2 for artificial viscosity calculations.

Consistency in the transport step is crucial for the conservation of mass and momentum. In nonuniform meshes, Equation (50) in Benítez-Llambay & Masset (2016) must be updated to

$$\begin{eqnarray}&&{v}_{y,{ijk}}^{n+1}=\displaystyle \frac{{{\rm{\Pi }}}_{y,{ijk}}^{-}{V}_{{ijk}}+{{\rm{\Pi }}}_{y,{ij}-1k}^{+}{V}_{{ij}-1k}}{{\rho }_{{ijk}}{V}_{{ijk}}+{\rho }_{{ij}-1k}{V}_{{ij}-1k}}.\end{eqnarray} \tag{ C5 }$$

The same must be done for v_x , v_z .