next previous
Up: FARGO: A fast eulerian


Subsections

3 New azimuthal transport algorithm

3.1 Overview

Let us take as an example the angular momentum conservation equation:

 \begin{displaymath}
\frac{\partial J}{\partial t}
+\underbrace{\frac{1}{r}\frac...
...}{\partial r}}_{\mbox{rad. transport}}
=\mbox{ Source terms }
\end{displaymath} (2)

where $J=\rho rv^\theta$. The transport equation of any HD quantity $\xi$ will look the same as the L.H.S. of Eq. (2).

Now without loss of generality we can rewrite Eq. (2) as:

 \begin{displaymath}
\frac{\partial J}{\partial t}
\!\!+\!\!\frac{1}{r}\frac{\pa...
...rtial(rv^{\rm r} J)}{\partial r}
\!=\!\! \mbox{ Source terms}
\end{displaymath} (3)

where u can be any quantity which does not depend on $\theta$. No assumption has been made on the behavior of J up to this point, and Eqs. (2) and (3) are strictly equivalent. If we take u to be the average azimuthal velocity $\overline v^\theta$, then Eq. (3) can be described as a composition of different steps, and each of them can be worked out independently with the well-known operator splitting technique:

A qualitative reason of why such a decomposition is valid is that the time evolution of the HD quantities can be described either by an observer sitting on a ring of radius r which rotates at any instant in time with the average azimuthal velocity, or by an observer at rest in an inertial frame. Now the time evolution of the system is of course observer-independent, which is why their observations are reconciled through the simple shift described by Eq. (5).

The idea on which the FARGO algorithm is based on is precisely to evolve the HD quantities through operators which mimic in a discrete way the different terms of Eq. (3). The source step, the radial transport step and the residual azimuthal velocity transport step are performed in a standard way (see e.g.Stone & Norman 1992). Now the last step in the operator-splitting described above, which corresponds to a simple shift which amounts to be $\overline v^\theta
\Delta t/r$ in one timestep, can be implemented in such a way that the matter can sweep an arbitrary number of cell widths in one timestep.

In order to lay down the basic mechanism by which FARGO works, let us take the following concrete example. We assume that, after the classical substeps (which are the source step, the radial transport and the residual azimuthal velocity transport), the material at a given radius r has to be shifted by 4.7 cells in one timestep (which means that $\overline v^\theta\Delta t/r\Delta\theta
=4.7$). What is actually done is that 4.7 is decomposed as 4.7=-0.3+5, i.e. the nearest integer and a remainder which by construction is lower or equal to 0.5 in absolute value. In the first substep of this shift step the material is shifted by this remainder (here -0.3), which can be achieved through a classical transport method since the remainder is lower or equal to 0.5 in absolute value (it has to be $\leq 1$ in order for the standard transport method to be possible), with the additional simplicity that the corresponding velocity field is uniform (which is actually why shift and transport happen to coincide in this special case, since there is no compression in the corresponding flow). The second substep just corresponds to an integer number of cells shift, which is done in our example simply by copying the content of cell j into cell j+5, for any j.

A more formal and detailed description of the FARGO algorithm is given in the next section.

3.2 Mathematical formulation of each step of the FARGO algorithm

In the modified algorithm, the azimuthal transport substep is split in several parts. We assume that the timestep $\Delta t$has already be chosen, and defer discussion of the timestep constraints until Sect. 3.3. We first compute the average azimuthal velocity at each radius:

 \begin{displaymath}
\overline{v}_i^\theta=\frac{1}{N_{\rm s}}\sum_{j=0}^{N_{\rm s}-1}
v_{ij}^{\theta a}.
\end{displaymath} (6)

We then introduce the residual velocity: $v_{ij}^{\theta {\rm res}}
=v_{ij}^{\theta a}-\overline{v}_i^\theta$, and the "shift number'' at each radius:

 \begin{displaymath}
n_i=E\left[\overline{v}^\theta_i\frac{\Delta t}{\Delta y_i}
\right]
\end{displaymath} (7)

where E[X] denotes the nearest integer to the real X. We define the constant residual velocity to be:

 \begin{displaymath}
v_i^{\theta{\rm cr}}=\overline{v}^\theta_i-n_i\frac{\Delta y_i}{\Delta t}\cdot
\end{displaymath} (8)

Hence the total velocity can be expressed as:

 \begin{displaymath}
v_{ij}^{\theta a}=v_i^{\theta {\rm SH}}+v_i^{\theta {\rm cr}}+
v_{ij}^{\theta {\rm res}}
\end{displaymath} (9)

where the "shift velocity'' $v_i^{\theta {\rm SH}}=
n_i\frac{\Delta y_i}{\Delta t}$ corresponds to a uniform shift of ni cells over one timestep.

We first transport the HD quantities according to the flow $v^{\theta {\rm res}}$:

 \begin{displaymath}
\xi_{ij}^c=\xi_{ij}^b+\frac{\Delta t}{\Delta y_i}
\left(\xi_...
...}^{b,*/v^{\theta {\rm res}}}v_{ij+1}^{\theta {\rm res}}\right)
\end{displaymath} (10)

then to the uniform flow $v^{\theta {\rm cr}}$:

 \begin{displaymath}
\xi_{ij}^d=\xi_{ij}^c+\frac{\Delta tv_i^{\theta {\rm cr}}}{\...
...theta {\rm cr}}}-
\xi_{ij+1}^{c,*/v^{\theta {\rm cr}}}\right).
\end{displaymath} (11)

We split the first part of the transport into two parts ( $v^{\theta {\rm res}}$ and $v^{\theta {\rm cr}}$) instead of using a single transport step with the velocity $v^{\theta {\rm res}}+v^{\theta {\rm cr}}$, in order to ensure (as can be checked below given the timestep constraints) that in each of these transport substeps the material sweeps at most half a cell (it could sweep up to one cell, but for reasons which will become clear in Sect. 4, we prefer to take a half cell limitation), and in order for the continuity considerations of Sect. 3.4 to apply. Finally, the quantities are transported along the $v^{\theta {\rm SH}}$uniform flow:

 \begin{displaymath}
\xi_{ij}^+=
\xi^d_{ij-n_i}.
\end{displaymath} (12)

Only the first two parts of this transport step introduce some numerical diffusion. The last one, given by Eq. (12), which in many cases corresponds to the largest part of the motion, does not introduce any numerical error, since it just corresponds to a circular permutation of the grid cells, or in other words it is just an integer discrete version of the shift given by Eq. (5).

A precise quantification of the lower numerical diffusivity of FARGO is beyond the scope of this paper though. An extremely rough estimation can be done in the case of the comparison of a standard method (in which the effective CFL ratio is a sizable fraction of one) and a FARGO method for which $n_i\ne 0$. If we assume that numerical effects will behave in azimuth as a physical viscosity would do, then the effective numerical viscosity in FARGO is about ni/C0 times lower than the standard method's one, where C0 is the CFL standard dimensionless limitation factor, which is detailed in the next section. Nevertheless a variety of numerical experiments can be found below which all show that FARGO's numerical diffusivity is smaller than the standard method's.

   
3.3 Timestep limitation

In the standard transport method, the timestep limitation arises from the combination of four different constraints (see e.g. Stone & Norman 1992), namely the fact that a flow advected test particle in cell [i,j] should not sweep a distance longer than $\Delta y_i$in azimuth nor longer than Ri+1-Ri in radius over one timestep (which introduces the limit timestep $\delta t_2$and $\delta t_3$ in Stone & Norman's paper), and that the wavefront of any wave present in the system should not travel across a whole cell over one timestep (Richtmyer & Morton 1957), which corresponds to the limit timestep $\delta t_1$ in Stone & Norman's paper. The last constraint comes from a stability limit arising from the viscosity (numerical or physical). With the modified azimuthal transport algorithm, the constraint on the azimuthal motion has to be modified slightly. Following Stone & Norman's notation, instead of writing $\delta t_3^{ij}=\Delta y_i/v^{\theta a}_{ij}$, we write:

 \begin{displaymath}
\delta t_{3}^{ij}
=\frac{\Delta y_i}
{v_{ij}^{\theta a}-\overline v_i^\theta}
=\frac{\Delta y_i}
{v_{ij}^{\theta \rm res}}
\end{displaymath} (13)

which means that the timestep limitation comes now from the perturbed azimuthal velocity, which results in a much higher absolute value of $\delta t_3$. Another limitation arises from the shear. Indeed we do not want the shear to disconnect the two neighboring cells [i,j] and [i+1,j] after one timestep. We write this condition as:

 \begin{displaymath}
\delta t_{\rm shear}^{ij}
={1 \over 2}{\left(
\frac{v_{ij}^{...
...}
-
\frac{v_{i+1j}^{\theta a}}{\Delta y_{i+1}}
\right)^{-1}
}.
\end{displaymath} (14)

Following Stone & Norman's notations, we finally adopt:
 
$\displaystyle \Delta t$ = $\displaystyle C_0/\left\{\max_{ij}\left[(\delta t_1^{ij})^{-2}+
(\delta t_2^{ij})^{-2}+
(\delta t_{3}^{ij})^{-2}\right.\right.$  
    $\displaystyle \left.\left. +(\delta t_4^{ij})^{-2}
+(\delta t_{\rm shear}^{ij})^{-2}\right]^{1/2}\right\}\cdot$ (15)

   
3.4 Continuity

At each timestep, $N_{\rm r}$ values of ni (with $i\in[0,N_{\rm r}-1]$), used in Eq. (12), are computed using Eq. (7). These integer values scale roughly as Ri-3/2. The shift on the central parts generally amounts to several cells over one timestep, while in the outer parts ni is small, and possibly zero. One can wonder whether or not problems may arise at the radii Ri where $n_i\neq n_{i-1}$ (i.e. at radii where the azimuthal shift corresponding to the third substep of the transport step is discontinuous). More generally we want to examine the question of the continuity of $\xi^+_{ij}$with respect to $\overline v_i^\theta\Delta t$. In order to check for this continuity, we assume $\overline v_i^\theta = \left(N+\frac{1}{2}
+\epsilon\right)
\frac{\Delta y_i}{\Delta t}$, where N is an integer, and we work out the behavior of $\xi_{ij}^+(\epsilon)$ in the vicinity of $\epsilon=0$. Since we have to use the explicit form of the " $*/v^{\theta a}$'' operator, we adopt the van Leer algorithm (van Leer 1977), which is widely used. Some straightforward algebra leads to:


\begin{displaymath}\xi_{ij}^+=\xi_{ij-N-1}^c+\left(\frac{1}{2}-\epsilon\right)
(\xi_{ij-N}^c-\xi_{ij-N-1}^c)\end{displaymath}


 \begin{displaymath}
\phantom{\xi_{ij}^+=}
-\left(\frac{1}{4}-\epsilon^2\right)
\frac{\Delta y_i}{2}({\rm d}\xi_{ij-N}^c-{\rm d}\xi_{ij-N-1}^c)
\end{displaymath} (16)

both for $\epsilon>0$ and $\epsilon<0$ provided $\vert\epsilon\vert < \frac{1}{2}$ and where the operator "d$\xi$'' is the van Leer slope. Equation (16) shows that the field $\xi_{ij}^+$is a continuous function of $\epsilon$ and hence of $\overline v_i^\theta$. In particular no special problem is to be expected from the discontinuities of niacross the disk.

3.5 Operators swapping

As we said in Sect. 2, it is a common practice to alternate the radial R and azimuthal T transport operators every other timestep. In this modified algorithm, Rshould usually be applied first, unless the velocity field is updated just after applying the T operator from the new momenta and new density fields, or unless special care is devoted to the j indices. Indeed swapping blindly the Rand T operators would result in moving radially the matter with the radial velocity it actually has $\sim n_i$ cells upwards, and would quickly end in a non-physical staggering everywhere $n_i\neq 0$.


next previous
Up: FARGO: A fast eulerian

Copyright The European Southern Observatory (ESO)