next previous
Up: FARGO: A fast eulerian


Subsections

   
4 Mono-dimensional tests

   
4.1 General considerations

In order to validate this modified transport algorithm, we present some 1D tests, and we compare the results of the standard method and of the FARGO method on a realistic test problem. We solve simultaneously the continuity and Navier Stokes equation for an isothermal gas (which has a non-vanishing but small kinematic viscosity):

 \begin{displaymath}
\frac{\partial\rho}{\partial t}+\frac{\partial(\rho v)}{\partial x}=0
\end{displaymath} (17)


 \begin{displaymath}
\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x...
...al \rho}{\partial x}+\nu
\frac{\partial^2v}{\partial x^2}\cdot
\end{displaymath} (18)

We assume that at rest the system has a uniform density $\rho_0$and sound speed $c_{\rm s}$. The waves which can propagate in this system have the following dispersion relationship:


\begin{displaymath}\omega=\pm\sqrt{k^2c_{\rm s}^2-\frac{k^4\nu^4}{4}}-i\frac{k^2\nu}{2}\end{displaymath}


 \begin{displaymath}
\mbox{or:~~~~}\omega=\pm kc_{\rm s}-i\frac{k^2\nu}{2}\mbox{~~~if~}\nu\ll\nu_{\rm lim}=\frac{2c_{\rm s}}{k}
\end{displaymath} (19)

which reduces to the standard dispersion relation for an undamped acoustic wave $\omega=\pm kc_{\rm s}$ provided the system is evolved for a time small compared to the damping timescale $\tau = \frac{2}{\nu k^2}$. This will be the case for the results we are going to present below, so that any apparent damping of the waves has a numerical origin. We do the following:
1.
We first analyze the propagation of a sound wave in the matter frame, i.e. we take as initial conditions:

 \begin{displaymath}
\rho(x)=s\rho_0\cos(kx)\mbox{~~~and~~~}v(x)=sc_{\rm s}\cos(kx)
\end{displaymath} (20)

where s is the wave relative amplitude. The polarization adopted corresponds to a rightwards propagating wave. According to Eq. (19), it propagates with a phase velocity which is $\Re \left(\frac{\omega}{k}\right)=c_{\rm s}$. We study this propagation with the standard transport algorithm (we are in the matter frame so there is no systematic average x-velocity, hence no need for a FARGO algorithm). We check that in this case the solution we get is accurate by varying the timestep and checking that the solution has converged.
2.
We then turn to a case where the setup is slightly modified. We take:

 \begin{displaymath}
\rho(x)=s\rho_0\cos(kx)\mbox{~~~and~~~}v(x)=v_0+sc_{\rm s}\cos(kx)
\end{displaymath} (21)

where v0 is a constant, which we choose much bigger than $c_{\rm s}$(which would correspond to the conditions of a thin keplerian disk, for example). The evolution of the system from this setup ought to be the same as before, since it merely corresponds to the same physical situation, but described from a frame moving at a constant speed -v0 wrt the first one, so one can invoke Galilean invariance to conclude that the wave profile evolution has to be the same. So any "good'' algorithm should approach as closely as possible the results of the matter frame simulations. We show that this is not quite the case with the standard transport method, which suffers from quite a high numerical dissipation, whereas FARGO behaves much better (not to mention its much faster execution). As a side result we also show that in this problem taking a CFL effective ratio (for the standard transport method) bigger than $\frac{1}{2}$ leads to an artificial and non-linear increase of the wave profile, and hence has to be avoided.

4.2 1D numerical results

We deal with a 1D grid composed of $N_{\rm s}=200$ cells, with periodic boundary conditions. The cell width is $\Delta x$ =0.0314, the isothermal sound speed is $c_{\rm s}=0.04$. The equilibrium density is $\Sigma_0=6\ 10^{-4}$. These parameters correspond roughly to the ones used in the numerical study of a protoplanet on a circular orbit at 5 AU embedded in a minimum mass protoplanetary disk (Hayashi et al.1985 or Bryden et al.1998), that are described in Sect. 5), when the central star mass and the protoplanet orbit radius are taken to be respectively the units of mass and distance. We present the results of different test runs in Fig. 1. The thick solid line represents the initial profile, which corresponds to a rightward propagating acoustic wave, with wavelength $\lambda=40\Delta x=1.256$. The relative amplitude of this sound wave is s=10-2. The thick dashed line represents the density profile at time t0=220, i.e. after the wave has traveled $c_{\rm s}t_0/\lambda=7$ times its own wavelength, when studied in the matter frame, i.e. when the velocity at t=0 is set to be only the perturbed velocity associated to the sound wave. The thick dashed profile is obtained with the standard transport algorithm (there is no need for the modified one in this case since we work in the matter frame), with a timestep $\Delta t=5\ 10^{-3}$. The curves obtained by choosing a much smaller timestep appear to coincide exactly with this one, hence we can consider this thick dashed line as the actual state the system must have at the date t0. This profile does not exactly coincide with the initial one because t0 is $\sim\frac{1}{7}$ of the profile steepening time $t_{\rm ps}\sim\frac{\lambda}{2c_{\rm s}s}$.


  \begin{figure}\par\psfig{file=zeroviscosity.ps,width=8.8cm}\par\end{figure} Figure 1: Compared evolution of an acoustic wave evolved with the standard transport algorithm and with the modified transport algorithm. We plot only two of the five wavelengths, i.e. 80 cells out of 200. Due to numerical effects the phase velocity of all these profiles do not exactly coincide with $c_{\rm s}$, so that after a time t0 their phases do not coincide. For this reason the profiles have been shifted so that they have all approximately the same phase in order to improve the clarity of the plot

Now if we just change the initial velocity by uniformly adding 1.0 to them at t=0, which means that we are no more in the matter frame, and we still work with the standard transport algorithm, then we get the dotted profile, which has $\sim 1/5$ the amplitude obtained from the computation in the matter frame. In this run the CFL ratio is $v\Delta t/\Delta x= 0.16$. In order to check the timestep dependency of this result, we redo this test with twice as smaller a timestep ( $\Delta t= 2.5\ 10^{-3}$) and we get the dash-dotted profile, which has about twice as smaller a density contrast than the previous curve. Note that if this effect were to be due to a physical kinematic viscosity $\nu$, then its value should be: $\nu\sim \frac{\lambda^2\log 5}{2\pi^2 t_0}\sim 5.8\ 10^{-4}$, much higher than the expected viscosity in a minimum mass protoplanetary disk ( $\nu \sim 10^{-5}$ in our dimensionless units). Now, instead of decreasing the timestep, we increase it and set $\Delta t=2.0\ 10^{-2}$ (hence the CFL ratio is about 0.64). We then get at time t0 the dot-dot-dot-dashed profile, which is not numerically damped but slightly amplified. With such a large timestep, we can use the modified transport algorithm, which in that case corresponds to a rightwards one cell shift and a leftwards normal transport with a remaining CFL ratio of 1-0.64=0.36. In that case we get the thin long-dashed profile. If we use the modified FARGO transport algorithm, we can still increase the timestep. The thin solid profile and the thin short-dashed profile have been obtained respectively with $\Delta t=4\ 10^{-2}$ (effective CFL ratio $\sim 1.3$) and $\Delta t=1.2\ 10^{-1}$ (effective CFL ratio $\sim 3.8$). We clearly see from these results that the FARGO transport algorithm leads to less numerical dissipation than the standard transport. From the first two tests in the non-comoving frame, one can conclude that increasing the number of timesteps over a given time interval with the standard transport algorithm increases the numerical dissipation (if the grid is moving wrt the matter frame with a velocity $v_0\neq 0$ and if the main part of the velocity comes from v0). A simple explanation for the lower numerical dissipation of the FARGO algorithm is that it requires less iterations as the timestep increases, and since most of the distance swept is achieved through an exact shift (a circular permutation), the numerical dissipation has to decrease as the timestep increases.


next previous
Up: FARGO: A fast eulerian

Copyright The European Southern Observatory (ESO)