Disco: A Moving Mesh Code for Gaseous Disks

Page Contents

  1. Code Overview
  2. Easy Hydro Tests
  3. Harder Hydro Tests
  4. Tests with Viscosity
  5. Planet-Disk Tests

‌

Code Overview

Disco solves hydrodynamic-like equations in conservation-law form:

$\partial_t u + \nabla \cdot \vec F = S$

It does so using a moving cylindrical mesh. The conserved quantities (u) are averaged over wedge-like volumes, and the fluxes (F) are averaged over space and time on the interfaces between adjacent volumes. Averaging the above conservation law over a wedge and over time gives:

$u^{n+1} dV^{n+1} = u^n dV^n - dt \Sigma_j F_j dA_j$

where the sum is over all faces attached to a given zone, and the index 'j' means the quantity is evaluated on face j. Now, if some of these faces are allowed to move with velocity $\vec w$, it is straightforward to show that this modifies the above equation in the following way:

$u^{n+1} dV^{n+1} = u^n dV^n - dt \Sigma_j ( F_j - u_j w \cdot \hat n ) dA_j$

where $\hat n$ is the unit normal to the face. The only remaining question is how to evaluate the $F_j's$ and $u_j's$, otherwise the entire method is essentially specified.

$\partial_t ( r \rho ) + \partial_r( r \rho v_r ) + (1/r)\partial_{\phi} ( r \rho v_{\phi} ) = 0$

$\partial_t ( r \rho v_r ) + \partial_r( r (\rho v_r^2 + P ) ) + (1/r)\partial_{\phi}( r \rho v_r v_{\phi} ) = \rho v_{\phi}^2 + P$

$\partial_t ( r^2 \rho v_{\phi} ) + \partial_r( r^2 \rho v_r v_{\phi} ) + (1/r)\partial_{\phi}( r^2 ( \rho v_{\phi}^2 + P ) ) = 0$

$\partial_t ( r ( \frac12 \rho \tilde v^2 + \epsilon ) ) \\ + \partial_r( r ( \frac12 \rho \tilde v^2 + \epsilon ) v_r + P v_r ) \\ + (1/r)\partial_{\phi}( r (\frac12 \rho \tilde v^2 + \epsilon) v_{\phi} + P \tilde v_{\phi}) \\ = r \rho v_r ( \Omega^2(r) r - \Omega'(r)(v_r/r - \Omega(r))r^2 )$

‌

Code Tests

  if( val < .1 ){
     bbb = 4.*(val+.125);
     ggg = 0.0; 
     rrr = 0.0; 
  }else if( val < .375){
     bbb = 1.0; 
     ggg = 4.*(val-.125);
     rrr = 0.0; 
  }else if( val < .625 ){
     bbb = 4.*(.625-val);
     rrr = 4.*(val-.375);
     ggg = bbb; 
     if( rrr > bbb ) ggg = rrr; 
  }else if( val < .875){
     bbb = 0.0; 
     ggg = 4.*(.875-val);
     rrr = 1.;
  }else{
     bbb = 0.0; 
     ggg = 0.0; 
     rrr = 4.*(1.125-val);
  }

‌

1D Tests

‌

Cylindrical Shock Tube

This test demonstrates Disco's ability to capture radially propagating shocks. For characteristics moving radially, Disco performs as a robust high-resolution shock-capturing code, even though the mesh motion is not employed.

It is also generally convenient to have a 2D shock solution to compare with numerically, which is why we have provided the numerical solution at high resolution.

Our domain extends from $0 < r < 1$.

We employ an adiabatic index of $\gamma = 5/3$.

Left State: $\rho = 1.0, P = 1.0$

Right State: $\rho = 0.125, P = 0.1$

We run this until a time $t = 0.25$.

Numerical Solution (8192 radial zones)

On the right we plot the density at this time, for a resolution of 128 radial zones.

‌

Cylindrical Isentropic Pulse

The isentropic pulse is a very simple nonlinear test for convergence of any code. We have set up a particular example and provided a high-resolution output which may be used as a basic code test. Convergence is easily checked by measuring entropy conservation. The initial setup is as follows:

$0 < r < 0.5$, uniform resolution.

$\rho = 1 + 3 e^{-80 r^2}$

$P = \rho^{\gamma}$

$v_r = v_\phi = 0$

For this problem, we choose an adiabatic index of $\gamma = 5/3$. The pulse explodes radially outward, and eventually forms a shock, but before the shock forms the equation $P = K \rho^{\gamma}$ continues to hold due to entropy conservation (this is straightforward to derive from Euler's equations; s = P/\rho^{\gamma} evolves as a conserved scalar, as long as the solution remains smooth).

We therefore calculate the error simply by verifying entropy conservation at a time $t = 0.1$, before the shock has formed:

$L_1 = { \int | P/\rho^{5/3} - 1.0 | dV \over \int dV}$.

We find fast convergence for this problem; for resolutions lower than 1024, we converge faster than second order. At higher resolutions, we converge at second order.

Numerical Solution (1024 radial zones)

L1 Error Data

This test is effectively 1D, as all motion is radial; it can be run with arbitrarily low angular resolution.

The test can also be used as a multidimensional test problem, simply by offsetting the origin of the pulse. We ran this problem with an origin offset by $\Delta y = 0.5$. These L1 errors are also included in the data file, though for these runs, the domain was extended out to $r = 1.0$, meaning the pulse is resolved at half the effective resolution of the centered case (also, the denominator in the definition of $L_1$ is quadrupled, which means the overall errors should be roughly the same in the offset case).

‌

2D Tests

‌

Smooth Vortex

This test helps to demonstrate convergence, and the effectiveness of the Disco code at preservation of contact discontinuities.

$0 < r < 5$, uniform resolution.

$~~\rho = 1.0$, $~~\Omega(r) = e^{-\frac12 r^2}$, $~~P(r) = 1 - \frac12 e^{-r^2}$.

The initial conditions are set so that centrifugal forces balance pressure gradients:

$\rho \Omega^2 r = \partial_r P$.

The vortex is trans-sonic (maximum Mach number of about 0.53). We calculate the error at $t = 10$ using the density:

$L_1 = { \int | \rho - 1.0 | dV \over \int dV}$.

L1 Error Data

Plotting this data, we see clear second-order convergence. In the plots on the right, we include a passive scalar to demonstrate the code's ability to maintain contact discontinuities and prevent artificial diffusion. When mesh motion is turned off, the contact discontinuity is significantly smeared out. With mesh motion turned on, we maintain the contact discontinuity very precisely. A calculation at a resolution of 256 radial zones is also included for comparison.

‌

Supersonic Keplerian Shear Flow

This is an extremely important test, as most of the real problems we wish to study have a Keplerian background flow. This stationary flow must be captured very accurately if we wish to study some subtle phenomenon which is added on top of this flow.

The setup for this problem is as follows:

$0.1 < r < 1.0$, Logarithmic radial zoning.

$~~\rho = 1.0$, $~~P = 0.01$, $~~\Omega(r) = 1/r^{3/2}$.

Boundary conditions are simply fixed at these initial conditions. A point mass is of course inserted at r = 0 so that centrifugal and gravitational forces are balanced. The equation of state is adiabatic, with index $\gamma = 5/3$. For this problem, we did not smooth the gravitational potential, since the point mass is not on the grid, and we require an exact solution (we could have instead modified the angular velocity to compensate, but this is simpler).

This gives a disk with a Mach number of 7.7 at the outer boundary, and 24.5 at the inner boundary.

Error is computed identically to the vortex problem. This is computed at a time $t = \pi$, after a half-orbit at the outer boundary, and about 16 orbits at the inner boundary. Truncation error generates transient waves which propagate radially, bouncing between the two boundaries. The L1 error is computed before these transient waves have fully dissipated. Convergence is very close to second-order. For 512 radial zones, the L1 error is at the $\text{few} \times 10^{-5}$ level.

L1 Error Data

Because this test is axisymmetric (and therefore effectively 1D), it is not very important that the zones move with the flow. However, for demonstration purposes we have added a passive scalar to the initial conditions:

$X = \theta(x)$ at $t = 0$

This passive scalar is plotted at time $t = 2\pi$, after the flow has had time to shear it into a spiral.

Because this problem is highly supersonic, most of the energy is kinetic, meaning that small relative errors in the energy density can lead to large errors in the pressure. For this reason, it is extremely beneficial to evolve a modified energy density:

$\frac12 \rho v^2 + \epsilon ~~ \rightarrow ~~ \frac12 \rho (\vec v - r \omega(r) \hat \phi)^2 + \epsilon$

where $\omega(r)$ is some prescribed analytic function. In our case, we set $\omega(r)$ to be the keplerian orbital velocity, so that we are mostly evolving thermal energy.

The tradeoff here is that energy is no longer conserved to machine precision, because the equations now have source terms involving $\omega(r)$ and its first derivative (which are known analytically). In practice, we do not conserve energy in Disco anyway, since gravity is introduced as a source term as well.

‌

Cylindrical Kelvin-Helmholtz Instability

Kelvin Helmholtz

$0 < r < 2$

if $r < 1$:

$~~~~\rho = 1$

$~~~~\Omega = 2$

$~~~~P = 10 + 2r^2$

if $r > 1$:

$~~~~\rho = 2$

$~~~~\Omega = 1$

$~~~~P = 11 + r^2$

$v_r = v_0 cos(6\phi) e^{-80(r-1)^2}$

$v_0 = 0.1$

$t = 2.0$

‌

Cartesian Shear Flow with Viscosity

In order to test Disco's implementation of viscosity, we need to perform a noncircular viscosity test. This is so that every term in the viscous equations is used. Fortunately, this is simple; all we need to do is set up a cartesian test problem and put it on our cylindrical grid. Of course, the cylindrical grid is certainly not ideal for this test, but the point of this test is to make sure all of our terms are implemented correctly, not to get extremely accurate results.

In cartesian coordinates, if we set up a flow with uniform density and pressure and with $\vec v = v(x) \hat y$, the Navier-Stokes equations reduce to:

$\dot v = \nu ~ v''$,

where $\nu$ is viscosity. This has the simple Green's function solution:

$v_y = { v_0 \over \sqrt{4 \pi \nu t}} exp( -{ (x-x_0)^2 \over 4 \nu t} ) $

We evolve this solution using the following parameters:

$0 < r < 2$

$~~\rho = 1$, $~~P = 1$, $~~x_0 = 1$, $~~v_0 = 0.001$, $~~\nu = 0.03$.

We start with the solution at time $t = 0.5$ and evolve until time $t = 1.0$. We used 64 radial zones, uniformly distributed. On the figure to the right we plot the solution at this time as a function of the x coordinate of each zone. Considering that the code is not designed for problems so misaligned with the grid, we capture the analytic solution reasonably well.

We also played with the various viscous fluxes and source terms in cylindrical coordinates, and found that all of these terms are necessary to capture this solution to this accuracy. Therefore, this test is an excellent way of being sure the viscosity is implemented correctly.

‌

Supersonic Keplerian Spreading

$0.2 < r < 1$

$\Omega(r) = 1/r^{3/2}$

$\Sigma(r) = 1 + 1 / \sqrt{r} + e^{-200(r-.5)^2}$

$P = 0.003$

$v_r = -\frac32{ \nu \over \Sigma r }$

$\nu = 10^{-5}$

$\gamma = 1.001$

This gives a disk orbiting at about Mach 45 in the vicinity of $r = 0.5$. We chose such a high Mach number because very thin disks are assumed when deriving the 1D diffusion equation for the surface density:

$\dot \Sigma = {3 \over r}( \sqrt{r} ( \sqrt{r} \Sigma \nu )' )'$

We wrote a simple 1D solver for this equation so that we can easily compare with Disco's performance on this problem. Of course, analytic Green's function solutions exist to this simple linear PDE, but I found it easier to solve the PDE numerically, that way I could easily choose whatever initial conditions I wanted and I didn't have to bother with any Bessel functions. For the initial conditions stated above for the surface density, we got the following solution at t = 100 by integrating this PDE:

1D Computed Solution

This is plotted on the panel to the right, compared with Disco's solution using 64 radial zones. The solutions do not seem to match identically, but this could be due to thin-disk assumptions made in deriving the 1D equation, or due to viscous heating, as we did not explicitly enforce isothermality. It could also possibly be due to the explicit form of the viscous terms in the hydro equations. At any rate, errors appear to be at the percent-level.

‌

Inviscid Disk with Earth-Mass Planet

It is important to be able to capture planet-disk interactions in a highly idealized context. We look at the linear interaction between an Earth-mass planet and a supersonic Keplerian disk whose Mach number is $\mathcal{M} = 20$ at the orbital radius. The initial conditions are not meant to mimic a protoplanetary disk, rather they are meant to produce an idealized environment which is easy for code comparison. We use a relatively small domain:

$0.5 < r < 1.5, ~~r_p = 1$

$\Omega(r) = 1/r^{3/2}, ~~\Sigma = 1, ~~P = 0.0025 = \rho/\mathcal{M}^2$

Here, to simplify the problem, we use a nearly isothermal equation of state:

$\gamma = 1.0001$

with boundary conditions fixed at the initial conditions. So far, this is just a supersonic Keplerian disk, similar to the test done above (except with a softer equation of state). The only additional ingredient is a point mass orbiting at $r_p = 1$.

We use the Earth-to-Sun mass ratio:

$q = 3 \times 10^{-6} = 0.024 q_{NL}$

where I have defined $q_{NL}$ according to the thermal mass threshold for nonlinearity $(q_{NL} = \mathcal{M}^{-3} = 1.25 \times 10^{-4})$.

The potential is of a point gravitating mass with smoothing length $\epsilon$:

$\phi(s) = { m_p \over \sqrt{s^2 + \epsilon^2} }$

where

$\epsilon = 0.025 = 0.5 h$

The simplest analytic comparison is given by the spiral wave produced by the planet. It is straightforward to calculate the shape of this wave from linear theory:

$\phi(r) = \phi_p + \text{sign}(r-r_p) (3 - 2\sqrt{r_p/r} - r/r_p) \mathcal{M}$

The spiral wave is quickly established in the first orbit. See figure.

‌

Something Else?

$\rho \rightarrow \rho ( 1 - (\gamma-1) \Phi / c^2 )^{1/(\gamma-1)}$

$P \rightarrow P ( 1 - (\gamma-1) \Phi / c^2 )^{\gamma/(\gamma-1)}$

‌