Page Contents

  1. Code Overview
  2. 1D Tests
  3. 2D Tests

Code Overview


Code Tests

1D Tests

Nonrelativistic Shock Tube

This is a very simple test problem which demonstrates the accurate shock-capturing of the code. This is important, since extremely strong shocks are expected in the code's main applications. The other nice thing about this test problem is it's a relatively simple set-up and therefore an easy first test for a code intended to solve Euler's equations in spherical coordinates. There is no analytic solution, so we provide a high-resolution result that can be used to test any code.

The set-up: $0 < r < 0.5$, uniform initial grid. Adiabatic index $\gamma = 5/3$.

if $r < 0.25$,

$~~\rho = P = 1$


$~~\rho = P = 0.1$.

We show the solution in the figure at time $t = 0.1$. An initial grid of 256 zones is compared with the high resolution result, and the shock and contact discontinuity are both captured with high accuracy. Note that the moving mesh causes zones to concentrate in the region between the shock and the contact discontinuity.

Hi-Res Solution (8192)

Isentropic Pulse

An isentropic wave can be used to easily test the convergence of any code. Any smooth isentropic initial conditions will remain isentropic as long as no shocks form. For this problem, we set up a spherical pulse which expands radially outward:

$0 < r < 0.5$

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

$P = \rho^{\gamma}$

In our case, we use the adiabatic index $\gamma = 5/3$. We can easily determine the convergence by calculating how precisely the code conserves entropy; that is, how closely the relationship $P = \rho^{\gamma}$ is maintained.

$L_1 = \int dr | P/\rho^{\gamma} - 1 |$

Hi-Res Solution (2048)

L1 Error Data

Sedov-Taylor Blastwave

$0 < r < 1.1$, uniform initial resolution

$\rho = 1$, $~~P = 10^{-6} + P_{exp} e^{-r^2/r_0^2}$

$r_0 = 0.02$, $~~P_{exp} = {\gamma - 1 \over \pi^{3/2} r_0^3 }$

For some reason we chose an adiabatic index $\gamma = 4/3$, even though this is a nonrelativistic problem. This gives an explosion with energy $E = 1.0$. By time $t = 1.0$, the state has evolved into the Sedov-Taylor solution. We found empirically that our solution matches closely with Sedov-Taylor as long as we displaced the solution a small amount in time to account for the finite radius of the initial conditions; we evaluate the analytic solution at time $t = 0.9872$. We have also included a code which evaluates this analytic solution (unfortunately, it is specific to the adiabatic index $\gamma = 4/3$ and I don't remember how I made the code in the first place so I didn't try to make it work for other values of $\gamma$).

Sedov Code

Adiabatic Expansion

We start with a gaussian bomb designed to have total energy and mass $E = M = 1.0$:

$\rho(r) = { e^{-r^2/R_0^2} \over \pi^{3/2} R_0^3 }$

$P(r) = (\gamma - 1) \rho(r)$

$R_0 = 10^{-8}$

We allow this solution to expand over about six orders of magnitude (until $t = 0.01$). At late times (but before the explosion has swept up a significant portion the ISM), the solution is described by a self-similar Hubble flow:

$v(r) = r/t$

$\rho(r) = {f(r/R) \over R^3}$

$P = \text{negligible}$.

where we have defined:

$R = v_{max} t$

$f(x)$ depends in the initial density distribution; we find for the nonrelativistic case that it simply emulates the initial condition, with some width we can determine empirically:

$f(x) = { e^{-x^2/a^2} \over \pi^{3/2} a^3 }$

Using conservation of energy we can calculate $v_{max}$:

$v_{max} = {2 \over \sqrt{3} a}$

We empirically determine for this problem

$a = 0.35 \rightarrow v_{max} = 3.3$

Spherical Collapse

$\dot r^2 = {GM \over r} - {GM \over r_0}$

${ dr \over \sqrt{ 1/r - 1/r_0 } } = - \sqrt{GM} dt $

$du = -{ dx \over \sqrt{ 1/x - 1 } } $

$x = r/r_0$

$u = t \sqrt{\frac83 \pi G \rho_0} $

$u = \sqrt{x(1-x)} - \frac12 \left( sin^{-1}( 2x-1 ) - \pi/2 \right)$

$t_{collapse} = {\pi/2 \over \sqrt{\frac83 \pi G \rho_0} }$

Blandford-McKee Blastwave

$E = 1.0, \rho_0 = 1.0$

$l = (E/\rho_0)^{1/3} = 1.0$

$\Gamma(t) = \sqrt{17 \over 8 \pi} (l / t)^{3/2}$

$\chi(r,t) = (1+8\Gamma^2)(1.-r/t)$

if $\chi > 1$ :

$~~\gamma(\chi) = \Gamma/\sqrt{2 \chi}$

$~~\rho(\chi) = 2 \rho_0 \Gamma^2 \chi^{-7/4}/\gamma(\chi)$

$~~P(\chi) = \frac23 \rho_0 \Gamma^2 \chi^{-17/12}$

Otherwise, density is set to $\rho_0$, velocity is set to zero and pressure is negligible. Note we have immediate simple scaling formulas for the maximum and average Lorentz factor:

$\gamma_{max}(t) = \Gamma/\sqrt{2} = \sqrt{17 \over 16 \pi} (l / t)^{3/2}$

$<\gamma \beta> = { \int r^2 dr (\gamma \beta) (\gamma^2 \rho h - P) \over \int r^2 dr (\gamma^2 \rho h - P) }$

$= { \int d\chi \gamma^3(\chi) P(\chi) \over \int d\chi \gamma^2(\chi) P(\chi) } = {17 \over 23}\gamma_{max}(t)$

We initiate the BMK solution at time $t_0 = .07$. At this time, the shock has Lorentz factor $\Gamma(t_0) = 44.4$. This of course is larger than the maximum Lorentz factor of the fluid, which is $\gamma_{max}(t_0) = 31.4$. The average Lorentz factor at this time is $<\gamma(t_0)> = 23.2$.

2D Tests

Nonrelativistic Rayleigh-Taylor Instability

$0.01 < r < 2.0$, initially logarithmic grid.

$0 < \theta < \pi/4$

$E = M = 1.0$

$\rho = {4 \pi/3 \over r_0^3}\theta(r_0-r) + 1$

$v = ( v_{max} r/r_0 ) \theta(r_0-r)$

$v_{max} = \sqrt{10/3}$

$r_0 = 0.1$

$t_0 = r_0 / v_{max} = 0.05477$

$t_0 < t < 2.0$