35.1 Hydrodynamics Test Problems

These problems are primarily designed to test the functioning of the hydrodynamics solvers within FLASH4.

35.1.1 Sod Shock-Tube

The Sod problem (Sod 1978) is a one-dimensional flow discontinuity problem that provides a good test of a compressible code's ability to capture shocks and contact discontinuities with a small number of cells and to produce the correct profile in a rarefaction. It also tests a code's ability to correctly satisfy the Rankine-Hugoniot shock jump conditions. When implemented at an angle to a multidimensional grid, it can be used to detect irregularities in planar discontinuities produced by grid geometry or operator splitting effects.

We construct the initial conditions for the Sod problem by establishing a planar interface at some angle to the $x$- and $y$-axes. The fluid is initially at rest on either side of the interface, and the density and pressure jumps are chosen so that all three types of nonlinear, hydrodynamic waves (shock, contact, and rarefaction) develop. To the “left” and “right” of the interface we have

$$\begin{array}{lclcclcl} \rho_{\rm L} &=& 1.0 & & & \rho_{\rm ... ...& 0.125 p_{\rm L} &=& 1.0 & & & p_{\rm R} &=& 0.1 \end{array}$$ (35.1)

The ratio of specific heats $\gamma$ is chosen to be 1.4 on both sides of the interface.

In FLASH, the Sod problem (Sod) uses the runtime parameters listed in Table 35.1 in addition to those supplied by default with the code. For this problem we use the Gamma equation of state alternative implementation and set gamma to 1.4. The default values listed in Table 35.1 are appropriate to a shock with normal parallel to the $x$-axis that initially intersects that axis at $x=0.5$ (halfway across a box with unit dimensions).

Table 35.1: Runtime parameters used with the Sod test problem.

Variable	Type	Default	Description
sim_rhoLeft	real	1	Initial density to the left of the interface ( $\rho_{\rm L}$)
sim_rhoRight	real	0.125	Initial density to the right ( $\rho_{\rm R}$)
sim_pLeft	real	1	Initial pressure to the left ($p_{\rm L}$)
sim_pRight	real	0.1	Initial pressure to the right ($p_{\rm R}$)
sim_uLeft	real	0	Initial velocity (perpendicular to interface) to the left ($u_{\rm L}$)
sim_uRight	real	0	Initial velocity (perpendicular to interface) to the right ($u_{\rm R}$)
sim_xangle	real	0	Angle made by interface normal with the $x$-axis (degrees)
sim_yangle	real	90	Angle made by interface normal with the $y$-axis (degrees)
sim_posn	real	0.5	Point of intersection between the interface plane and the $x$-axis

Figure: Comparison of numerical and analytical solutions to the Sod problem. A 2D grid with six levels of refinement is used. The shock normal is parallel to the $x$-axis.

Figure: Comparison of numerical solutions to the Sod problem for two different angles ($\theta$) of the shock normal relative to the $x$-axis. A 2D grid with six levels of refinement is used.

Figure 35.2 shows the result of running the Sod problem with FLASH on a two-dimensional grid with the analytical solution shown for comparison. The hydrodynamical algorithm used here is the directionally split piecewise-parabolic method (PPM) included with FLASH. In this run the shock normal is chosen to be parallel to the $x$-axis. With six levels of refinement, the effective grid size at the finest level is $256^2$, so the finest cells have width 0.00390625. At $t=0.2$, three different nonlinear waves are present: a rarefaction between $x = 0.263$ and $x = 0.486$, a contact discontinuity at $x = 0.685$, and a shock at $x = 0.850$. The two discontinuities are resolved with approximately two to three cells each at the highest level of refinement, demonstrating the ability of PPM to handle sharp flow features well. Near the contact discontinuity and in the rarefaction, we find small errors of about $1-2\%$ in the density and specific internal energy, with similar errors in the velocity inside the rarefaction. Elsewhere, the numerical solution is close to exact; no oscillations are present.

Figure 35.3 shows the result of running the Sod problem on the same two-dimensional grid with different shock normals: parallel to the $x$-axis ( $\theta=0^\circ$) and along the box diagonal ( $\theta=45^\circ$). For the diagonal solution, we have interpolated values of density, specific internal energy, and velocity to a set of 256 points spaced exactly as in the $x$-axis solution. This comparison shows the effects of the second-order directional splitting used with FLASH on the resolution of shocks. At the right side of the rarefaction and at the contact discontinuity, the diagonal solution undergoes slightly larger oscillations (on the order of a few percent) than the $x$-axis solution. Also, the value of each variable inside the discontinuity regions differs between the two solutions by up to 10%. However, the location and thickness of the discontinuities is the same for the two solutions. In general, shocks at an angle to the grid are resolved with approximately the same number of cells as shocks parallel to a coordinate axis.

Figure 35.4 presents a colormap plot of the density at $t=0.2$ in the diagonal solution together with the block structure of the AMR grid. Note that regions surrounding the discontinuities are maximally refined, while behind the shock and contact discontinuity, the grid has de-refined, because the second derivative of the density has decreased in magnitude. Because zero-gradient outflow boundaries were used for this test, some reflections are present at the upper left and lower right corners, but at $t=0.2$ these have not yet propagated to the center of the grid.

Figure: Density in the diagonal 2D Sod problem with six levels of refinement at $t=0.2$. The outlines of AMR blocks are shown (each block contains $8\times 8$ cells).

SodStep Example: FLASH4 also contains under the name SodStep a variant of the Sod problem. This setup is provided as an example for setting up simulations on domains with steps and obstacles. See the files in the SodStep simulation directory and the Simulation_defineDomain description for more information on how to use this feature.

35.1.2 Variants of the Sod Problem in Curvilinear Geometries

Variants of the Sod problems can be set up in in various geometries in order to test the handling of non-Cartesion geometries.

• An axisymmetric variant of the Sod problem can be configured by setting up the regular Sod simulation with ./setup Sod -auto -2d -geometry=cylindrical and using runtime parameters that include geometry = "cylindrical". Use sim_xangle = 0 to configure an initial shock front that is shaped like a cylinder. Results as in those discussed in Toro 1999 can be obtained.
• A spherically symmetric variant of the Sod problem can be configured by setting up the regular Sod simulation with ./setup Sod -auto -1d -geometry=spherical and using runtime parameters that include geometry = "spherical". Again results as in those discussed in Toro 1999 can be obtained.
• To test the behavior of FLASH solutions when the physical symmetry of the problem does not match the geometry of the simulation, a separate simulation is provided under the name SodSpherical. To use this, configure with ./setup SodSpherical -auto -2d -geometry=spherical and using runtime parameters that include geometry = "spherical". As a 2D setup, SodSpherical represents physically axisymmetric initial conditions in spherical coordinates. The physical problem can be chosen to be the same as in the previous case with cylindrical Sod. Again results as in those discussed in Toro 1999 can be obtained.
• The SodSpherical setup can also configured in 1D and will act like the 1D Sod setup in that case.

35.1.3 Interacting Blast-Wave Blast2

This Blast2 problem was originally used by Woodward and Colella (1984) to compare the performance of several different hydrodynamical methods on problems involving strong shocks and narrow features. It has no analytical solution (except at very early times), but since it is one-dimensional, it is easy to produce a converged solution by running the code with a very large number of cells, permitting an estimate of the self-convergence rate. For FLASH, it also provides a good test of the adaptive mesh refinement scheme.

The initial conditions consist of two parallel, planar flow discontinuities. Reflecting boundary conditions are used. The density is unity and the velocity is zero everywhere. The pressure is large at the left and right and small in the center

$$\begin{array}{lclcclclcclcl} p_{\rm L} &=& 1000, && p_{\rm M} &=& 0.01, && p_{\rm R} &=& 100 . \end{array}$$ (35.2)

The equation of state is that of a perfect gas with $\gamma=1.4$.

Figure 35.5 shows the density and velocity profiles at several different times in the converged solution, demonstrating the complexity inherent in this problem. The initial pressure discontinuities drive shocks into the middle part of the grid; behind them, rarefactions form and propagate toward the outer boundaries, where they are reflected back into the grid. By the time the shocks collide at $t=0.028$, the reflected rarefactions have caught up to them, weakening them and making their post-shock structure more complex. Because the right-hand shock is initially weaker, the rarefaction on that side reflects from the wall later, so the resulting shock structures going into the collision from the left and right are quite different. Behind each shock is a contact discontinuity left over from the initial conditions (at $x\approx 0.50$ and 0.73). The shock collision produces an extremely high and narrow density peak. The peak density should be slightly less than 30. However, the peak density shown in Figure 35.5 is only about 18, since the maximum value of the density does not occur at exactly $t=0.028$. Reflected shocks travel back into the colliding material, leaving a complex series of contact discontinuities and rarefactions between them. A new contact discontinuity has formed at the point of the collision ( $x\approx 0.69$). By $t=0.032$, the right-hand reflected shock has met the original right-hand contact discontinuity, producing a strong rarefaction, which meets the central contact discontinuity at $t=0.034$. Between $t=0.034$ and $t=0.038$, the slope of the density behind the left-hand shock changes as the shock moves into a region of constant entropy near the left-hand contact discontinuity.

Figure 35.5: Density and velocity profiles in the Woodward-Colella interacting blast-wave problem Blast2 as computed by FLASH using ten levels of refinement.

Figure 35.6 shows the self-convergence of density and pressure when FLASH is run on this problem. We compare the density, pressure, and total specific energy at $t=0.038$ obtained using FLASH with ten levels of refinement to solutions using several different maximum refinement levels. This figure plots the L1 error norm for each variable $u$, defined using

$${\cal E}(N_{\rm ref};u) \equiv {1\over N(N_{\rm ref})} \sum_{i=1}... ..._{\rm ref})}\left\vert{u_i(N_{\rm ref}) - A u_i(10)\over u_i(10)}\right\vert ,$$ (35.3)

against the effective number of cells ( $N(N_{\rm ref})$). In computing this norm, both the `converged' solution $u(10)$ and the test solution $u(N_{\rm ref})$ are interpolated onto a uniform mesh having $N(N_{\rm ref})$ cells. Values of $N_{\rm ref}$ between 2 (corresponding to cell size $\Delta x=1/16$) and 9 ( $\Delta x=1/2048$) are shown.

Although PPM is formally a second-order method, the convergence rate is only linear. This is not surprising, since the order of accuracy of a method applies only to smooth flow and not to flows containing discontinuities. In fact, all shock capturing schemes are only first-order accurate in the vicinity of discontinuities. Indeed, in their comparison of the performance of seven nominally second-order hydrodynamic methods on this problem, Woodward and Colella found that only PPM achieved even linear convergence; the other methods were worse. The error norm is very sensitive to the correct position and shape of the strong, narrow shocks generated in this problem.

The additional runtime parameters supplied with the 2blast problem are listed in Table 35.2. This problem is configured to use the perfect-gas equation of state (gamma) with gamma set to 1.4 and is run in a two-dimensional unit box. Boundary conditions in the $y$-direction (transverse to the shock normals) are taken to be periodic.

Table 35.2: Runtime parameters used with the 2blast test problem.

Variable	Type	Default	Description
rho_left	real	1	Initial density to the left of the left interface ( $\rho_{\rm L}$)
rho_mid	real	1	Initial density between the two interfaces ( $\rho_{\rm M}$)
rho_right	real	1	Initial density to the right of the right interface ( $\rho_{\rm R}$)
p_left	real	1000	Initial pressure to the left ($p_{\rm L}$)
p_mid	real	0.01	Initial pressure in the middle ($p_{\rm M}$)
p_right	real	100	Initial pressure to the right ($p_{\rm R}$)
u_left	real	0	Initial velocity (perpendicular to interface) to the left ($u_{\rm L}$)
u_mid	real	0	Initial velocity (perpendicular to interface) in the middle ($u_{\rm M}$)
u_right	real	0	Initial velocity (perpendicular to interface) to the right ($u_{\rm R}$)
xangle	real	0	Angle made by interface normal with the $x$-axis (degrees)
yangle	real	90	Angle made by interface normal with the $y$-axis (degrees)
posnL	real	0.1	Point of intersection between the left interface plane and the $x$-axis
posnR	real	0.9	Point of intersection between the right interface plane and the $x$-axis

Figure 35.6: Self-convergence of the density, pressure, and total specific energy in the Blast2 test problem.

35.1.4 Sedov Explosion

The Sedov explosion problem (Sedov 1959) is another purely hydrodynamical test in which we check the code's ability to deal with strong shocks and non-planar symmetry. The problem involves the self-similar evolution of a cylindrical or spherical blast wave from a delta-function initial pressure perturbation in an otherwise homogeneous medium. We provide two different ways to generate the initial conditions:

1. We deposit a quantity of energy $E=1$ into a small region of radius $\delta r$ at the center of the grid. The pressure inside this volume $p_0'$ is given by

   $$p_0' = {3(\gamma-1)E\over(\nu+1)\pi \delta r^\nu} ,$$ (35.4)

   
   
   where $\nu=2$ for cylindrical geometry and $\nu=3$ for spherical geometry.
2. We initialize from a precomputed pseudo-analytical solution by reading in a 1D profile.

We set the ratio of specific heats $\gamma=1.4$. In running this problem we choose $\delta r$ to be 3.5 times as large as the finest adaptive mesh resolution in order to minimize effects due to the Cartesian geometry of our grid. The density is set equal to $\rho_0=1$ everywhere, and the pressure is set to a small value $p_0=10^{-5}$ everywhere but in the center of the grid. The fluid is initially at rest. In the self-similar blast wave that develops for $t>0$, the density, pressure, and radial velocity are all functions of $\xi \equiv r/R(t)$, where

$$R(t) = C_\nu(\gamma) \left({Et^2\over \rho_0}\right)^{1/(\nu+2)} .$$ (35.5)

Here $C_\nu$ is a dimensionless constant depending only on $\nu$ and $\gamma$; for $\gamma=1.4$, $C_2 \approx C_3 \approx 1$ to within a few percent. Just behind the shock front at $\xi = 1$ the analytical solution is

$$\rho = \rho_1 \equiv {\gamma+1\over\gamma-1}\rho_0$$

$$p = p_1 \equiv {2\over\gamma+1}\rho_0 u^2$$ (35.6)

$$v = v_1 \equiv {2\over\gamma+1}u ,$$

where $u \equiv dR/dt$ is the speed of the shock wave. Near the center of the grid,

$$\rho(\xi)/\rho_1 \propto \xi^{\nu/(\gamma-1)}$$

$$p(\xi)/p_1 = {\rm constant}$$ (35.7)

$$v(\xi)/v_1 \propto \xi .$$

Figure: Comparison of numerical and analytical solutions to the Sedov problem in two dimensions. Numerical solution values are averages in radial bins at the finest AMR grid resolution $N_{\rm ref}$ in each run.

Figure 35.7 shows density, pressure, and velocity profiles in the two-dimensional, cylindrical Sedov problem at $t=0.05$. Solutions obtained with FLASH on grids with 2, 4, 6, and 8 levels of refinement are shown in comparison with the analytical solution. In this figure we have computed average radial profiles in the following way. We interpolated solution values from the adaptively gridded mesh used by FLASH onto a uniform mesh having the same resolution as the finest AMR blocks in each run. Then, using radial bins with the same width as the cells in the uniform mesh, we binned the interpolated solution values, computing the average value in each bin. At low resolutions, errors show up as density and velocity overestimates behind the shock, underestimates of each variable within the shock, and a very broad shock structure. However, the central pressure is accurately determined, even for two levels of refinement. Because the density goes to a finite value rather than to its correct limit of zero, this corresponds to a finite truncation of the temperature (which should go to infinity as $r\rightarrow 0$). This error results from depositing the initial energy into a finite-width region rather than starting from a delta function. As the resolution improves and the value of $\delta r$ decreases, the artificial finite density limit also decreases; by $N_{\rm ref}=6$ it is less than 0.2% of the peak density. Except for the $N_{\rm ref}=2$ case, which does not show a well-defined peak in any variable, the shock itself is always captured with about two cells. The region behind the shock containing 90% of the swept-up material is represented by four cells in the $N_{\rm ref}=4$ case, 17 cells in the $N_{\rm ref}=6$ case, and 69 cells for $N_{\rm ref}=8$. However, because the solution is self-similar, for any given maximum refinement level, the structure will be four cells wide at a sufficiently early time. The behavior when the shock is under-resolved is to underestimate the peak value of each variable, particularly the density and pressure.

Figure 35.8 shows the pressure field in the 8-level calculation at $t=0.05$ together with the block refinement pattern. Note that a relatively small fraction of the grid is maximally refined in this problem. Although the pressure gradient at the center of the grid is small, this region is refined because of the large temperature gradient there. This illustrates the ability of PARAMESH to refine grids using several different variables at once.

Figure: Pressure field in the 2D Sedov explosion problem with 8 levels of refinement at $t=0.05$. The outlines of the AMR blocks are overlaid on the pressure colormap.

We have also run FLASH on the spherically symmetric Sedov problem in order to verify the code's performance in three dimensions. The results at $t=0.05$ using five levels of grid refinement are shown in Figure 35.9. In this figure we have plotted the average values as well as the root-mean-square (RMS) deviations from the averages. As in the two-dimensional runs, the shock is spread over about two cells at the finest AMR resolution in this run. The width of the pressure peak in the analytical solution is about 1 1/2 cells at this time, so the maximum pressure is not captured in the numerical solution. Behind the shock the numerical solution average tracks the analytical solution quite well, although the Cartesian grid geometry produces RMS deviations of up to 40% in the density and velocity in the de-refined region well behind the shock. This behavior is similar to that exhibited in the two-dimensional problem at comparable resolution.

Figure: Comparison of numerical and analytical solutions versus radius $r$ to the spherically symmetric Sedov problem. A 3D grid with five levels of refinement is used.

The additional runtime parameters supplied with the Sedov problem are listed in Table 35.3. This problem is configured to use the perfect-gas equation of state (Gamma) with gamma set to 1.4. It is simulated in a unit-sized box.

Table 35.3: Runtime parameters used with the Sedov test problem.

Variable	Type	Default	Description
sim_pAmbient	real	$10^{-5}$	Initial ambient pressure ($p_0$)
sim_rhoAmbient	real	1	Initial ambient density ($\rho_0$)
sim_expEnergy	real	1	Explosion energy ($E$)
sim_rInit	real	0.05	Radius of initial pressure perturbation ($\delta r$)
sim_xctr	real	0.5	$x$-coordinate of explosion center
sim_yctr	real	0.5	$y$-coordinate of explosion center
sim_zctr	real	0.5	$z$-coordinate of explosion center
sim_nSubZones	integer	7	Number of sub-cells in cells for applying the 1D profile

35.1.4.1 Sedov Self-Gravity

Another variant of the Sedov problem is included in the release which runs with spherical coordinates in one dimension. The Sedov Self-Gravity problem also allows the effects of gravitational acceleration where the gravitational potential is calculated using the multipole solver. Figure 35.10 and 35.11 show the snapshots of the density profile and the gravitational potential at two different times during the evolution. The first snapshot is at $t=0.5$ sec, when evolution is halfway through, while the second snapshot is at the end of the evolution, where $t=1.0$ sec.

Figure 35.10: Snapshots of Sedov Self-gravity density profile and gravitational potential at time t=0.5 sec.

Figure 35.11: Snapshots of Sedov Self-gravity density profile and gravitational potential at time t=1.0 sec.

35.1.5 Isentropic Vortex

The two-dimensional isentropic vortex problem is often used as a benchmark for comparing numerical methods for fluid dynamics. The flow-field is smooth (there are no shocks or contact discontinuities) and contains no steep gradients, and the exact solution is known. It was studied by Yee, Vinokur, and Djomehri (2000) and by Shu (1998). In this subsection the problem is described, the FLASH control parameters are explained, and some results demonstrating how the problem can be used are presented.

The simulation domain is a square, and the center of the vortex is located at $(x_{ctr}, y_{ctr})$. The flow-field is defined in coordinates centered on the vortex center $(x' = x - x_{ctr}, y' = y - y_{ctr})$ with $r^2 = {x'}^2 + {y'}^2$. The domain is periodic, but it is assumed that off-domain vortexes do not interact with the primary; practically, this assumption can be satisfied by ensuring that the simulation domain is large enough for a particular vortex strength. We find that a domain size of $10 \times 10$ (specified through the Grid runtime parameters xmin, xmax, ymin, and ymax) is sufficiently large for a vortex strength (defined below) of 5.0. In the initialization below, $x'$ and $y'$ are the coordinates with respect to the nearest vortex in the periodic sense.

The ambient conditions are given by $\rho_\infty$, $u_\infty$, $v_\infty$, and $p_\infty$, and the non-dimensional ambient temperature is $T^*_\infty = 1.0$. Using the equation of state, the (dimensional) $T_\infty$ is computed from $p_\infty$ and $\rho_\infty$. Perturbations are added to the velocity and nondimensionalized temperature, $u = u_\infty + \delta u$, $v = v_\infty + \delta v$, and $T^* = T^*_\infty + \delta T^*$ according to

$$\delta u = -y' {\frac {\beta} {2 \pi}} \exp \left( {\frac {1-r^2} {2}} \right)$$ (35.8)

$$\delta v = x' {\frac {\beta} {2 \pi}} \exp \left( {\frac {1-r^2} {2}} \right)$$ (35.9)

$$\delta T^* = - { \frac {(\gamma - 1 ) \beta} {8 \gamma \pi^2}} \exp \left( {1-r^2} \right) ,$$ (35.10)

where $\gamma=1.4$ is the ratio of specific heats and $\beta=5.0$ is a measure of the vortex strength. The temperature and density are then given by

$$T = {\frac{T_\infty}{T^*_\infty} } T^*$$ (35.11)

$$\rho = \rho_\infty \left( {\frac{T}{T_\infty} } \right)^{\frac{1}{\gamma-1} } .$$ (35.12)

At any location in space, the conserved variables (density, $x$- and $y$-momentum, and total energy) can be computed from the above quantities. The flow-field is initialized by computing cell averages of the conserved variables; in each cell, the average is approximated by averaging over ${\tt nx\_subint} \times {\tt ny\_subint}$ subintervals. The runtime parameters for the isentropic vortex problem are listed in Table 35.4.

Table 35.4: Parameters for the IsentropicVortex problem.

Variable	Type	Default	Description
p_ambient	real	1.0	Initial ambient pressure ( $p_{\infty}$)
rho_ambient	real	1.0	Initial ambient density ( $\rho_{\infty}$)
u_ambient	real	1.0	Initial ambient $x$-velocity ( $u_{\infty}$)
v_ambient	real	1.0	Initial ambient $y$-velocity ( $v_{\infty}$)
vortex_strength	real	5.0	Non-dimensional vortex strength
xctr	real	0.0	$x$-coordinate of vortex center
yctr	real	0.0	$y$-coordinate of vortex center
nx_subint	integer	10	number of subintervals in $x$-direction
ny_subint	integer	10	number of subintervals in $y$-direction

Figure 35.12 shows the exact density distribution represented on a $40 \times 40$ uniform grid with $-5.0 \leq x, y \leq 5.0$. The borders of each grid block ( $8\times 8$ cells) are superimposed. In addition to the shaded representation, contour lines are shown for $\rho = 0.95$, 0.85, 0.75, and 0.65. The density distribution is radially symmetric, and the minimum density is $\rho_{min} = 0.510287$. Because the exact solution of the isentropic vortex problem is the initial solution shifted by $(u_\infty t, v_\infty t)$, numerical phase (dispersion) and amplitude (dissipation) errors are easy to identify. Dispersive errors distort the shape of the vortex, breaking its symmetry. Dissipative errors smooth the solution and flatten extrema; for the vortex, the minimum in density at the vortex core will increase.

Figure: Density at $t=0.0$ for the isentropic vortex problem. Shown are the initial condition and the exact solution at $t=10.0, 20.0, \ldots$.

A numerical simulation using the PPM scheme was run to illustrate such errors. The simulation used the same grid shown in Figure 35.12 with the same contour levels and color values. The grid is intentionally coarse and the evolution time long to make numerical errors visible. The vortex is represented by approximately 8 grid points in each coordinate direction and is advected diagonally with respect to the grid. At solution times of $t=10, 20, \ldots$, etc., the vortex should be back at its initial location.

Figure 35.13 shows the solution at $t=50.0$; only slight differences are observed. The density distribution is almost radially symmetric, although the minimum density has risen to $0.0537360$. Accumulating dispersion error is clearly visible at $t=100.0$ (Figure 35.14), and the minimum density is now $0.601786$.

Figure: Density at $t=50.0$ for the isentropic vortex problem.

Figure: Density at $t=100.0$ for the isentropic vortex problem.

Figure 35.15 shows the density near $y=0.0$ at three simulation times. The black line shows the initial condition. The red line corresponds to $t=50.0$ and the blue line to $t=100.0$. For the later two times, the density is not radially symmetric. The lines plotted are just representative profiles for those times, which give an idea of the magnitude and character of the errors.

Figure: Representative density profiles for the isentropic vortex near $y=0.0$ at $t=0.0$ (black), $t=50.0$ (red), and $t=100.0$ (blue).

35.1.6 The double Mach reflection problem

This numerical planar shock test problem introduced by Woodward and Colella (1984) simulates an evolution of an unsteady planar shock that is incident on an oblique surface. Initially, the incident planar shock begins to propagate to the bottom surface at a $30^{\circ}$ angle with the shock Mach number of 10. Later, the solution to this problem produces a self-similar wave pattern that corresponds to the double Mach reflection. Following many other numerical setups of this problem, we tilt the incident shock rather than the reflecting wall so as to avoid numerical complications in modeling an oblique physical boundary.

The initial setup involves a Mach 10 shock in air, $\gamma=1.4$, on a rectangular 2D domain of size $[0,4]\times[0,1]$. The reflecting wall is represented as the bottom surface of the domain, beginning at $x=1/6$. The velocity normal to the incident shock in the post-shock region is 8.25, and the flow is at rest in the pre-shock region. The undisturbed air ahead of the shock has a density of 1.4 and a pressure of 1. See the initial density profile in Figure 35.16 resolved on 6 refinement AMR levels using 16$^2$-cell block size.

The boundary condition on $0\le x \le 1/6$ at $y=0$ is fixed in time with the initial values so that the reflected shock is attached to the bottom surface. We impose reflecting boundary condition (i.e., negating the $y$ velocity field $v$) on the rest of the bottom surface. On the top surface of $y=1$, we allow the initial Mach 10 shock proceeds exactly as a function of time in order that the numerical evolution follows the oblique shock propagation without any planar distortion. At $x = 0$, we impose a supersonic inflow boundary condition, and the outflow condition at $x=4$.

In later time, the solution develops to form two Mach stems and two contact discontinuities, as shown in Figure 35.17 the density at $t=2.5$ sec. Also shown in Figure 35.18 are 30 levels of contour lines of pressure, illustrating the evolved flow discontinuities at $t=2.5$ sec. We can also see that the numerical solution is smooth and non-oscillatory in the region bounded by the curved reflected shock, the bottom surface and the second Mach stem (see more discussion in Woodward and Colella, 1984).

The 5th order WENO method in the unsplit hydrodynamics solver is used with a choice of hybrid Riemann solver which selectively adopts HLL only at strong shocks and HLLC otherwise.

Figure: The initial density at $t=0$ visualized with 6 levels of AMR block structures.

Figure: Density profile at $t=2.5$. Two contact discontinuities are denoted as "CD", along with two Mach stems, a curved reflected shock, and a jet formulation of the denser fluid along the bottom surface.

Figure: The 30 levels of contour lines of pressure at $t=2.5$.

Table 35.5: Runtime parameters used with the double Mach reflection test problem.

Variable	Type	Default	Description
sim_rhoLeft	real	8	Initial density to the left of the shock ( $\rho_{\rm L}$)
sim_rhoRight	real	1.4	Initial density to the right ( $\rho_{\rm R}$)
sim_pLeft	real	116.5	Initial pressure to the left ($p_{\rm L}$)
sim_pRight	real	1	Initial pressure to the right ($p_{\rm R}$)
sim_uLeft	real	7.1447096	Initial $x$-velocity to the left ($u_{\rm L}$)
sim_uRight	real	0	Initial $x$-velocity to the right ($u_{\rm R}$)
sim_vLeft	real	-4.125	Initial $y$-velocity to the left ($v_{\rm L}$)
sim_vRight	real	0	Initial $y$-velocity to the right ($v_{\rm R}$)
sim_xangle	real	-30	Angle made by shock normal with the $x$-axis (degrees)
sim_yangle	real	90	Angle made by shock normal with the $y$-axis (degrees)
sim_posn	real	1/6	Point of intersection between the shock plane and the $x$-axis

35.1.7 Wind Tunnel With a Step

The problem of a wind tunnel containing a step, WindTunnel was first described by Emery (1968), who used it to compare several hydrodynamical methods. Woodward and Colella (1984) later used it to compare several more advanced methods, including PPM. Although it has no analytical solution, this problem is useful because it exercises a code's ability to handle unsteady shock interactions in multiple dimensions. It also provides an example of the use of FLASH to solve problems with irregular boundaries.

The problem uses a two-dimensional rectangular domain three units wide and one unit high. Between $x=0.6$ and $x=3$ along the $x$-axis is a step 0.2 units high. The step is treated as a reflecting boundary, as are the lower and upper boundaries in the $y$-direction. For the right-hand $x$-boundary, we use an outflow (zero gradient) boundary condition, while on the left-hand side we use an inflow boundary. In the inflow boundary cells, we set the density to $\rho_0$, the pressure to $p_0$, and the velocity to $u_0$, with the latter directed parallel to the $x$-axis. The domain itself is also initialized with these values. We use

$$\rho_0 = 1.4, \qquad p_0 = 1, \qquad u_0 = 3 , \qquad \gamma = 1.4,$$ (35.13)

which corresponds to a Mach 3 flow. Because the outflow is supersonic throughout the calculation, we do not expect reflections from the right-hand boundary.

Figure 35.19: Density and velocity in the Emery wind tunnel test problem, as computed with FLASH. A 2D grid with five levels of refinement is used.

Figure 35.19: Density and velocity in the Emery wind tunnel test problem (continued).

The additional runtime parameters supplied with the WindTunnel problem are listed in Table 35.6. This problem is configured to use the perfect-gas equation of state (Gamma) with gamma set to 1.4. We also set xmax$=3$, ymax$=1$, Nblockx$=15$, and Nblocky$=5$ in order to create a grid with the correct dimensions. The version of Simulation_defineDomain supplied with this problem removes all but the first three top-level blocks along the lower edge of the grid to generate the step, and gives REFLECTING boundaries to the obstacle blocks. Finally, we use xl_boundary_type = "user" (USER_DEFINED condition) and xr_boundary_type = "outflow" (OUTFLOW boundary) to instruct FLASH to use the correct boundary conditions in the $x$-direction. Boundaries in the $y$-direction are reflecting (REFLECTING).

Until $t=12$, the flow is unsteady, exhibiting multiple shock reflections and interactions between different types of discontinuities. Figure 35.19 shows the evolution of density and velocity between $t=0$ and $t=4$ (the period considered by Woodward and Colella). Immediately, a shock forms directly in front of the step and begins to move slowly away from it. Simultaneously, the shock curves around the corner of the step, extending farther downstream and growing in size until it strikes the upper boundary just after $t=0.5$. The corner of the step becomes a singular point, with a rarefaction fan connecting the still gas just above the step to the shocked gas in front of it. Entropy errors generated in the vicinity of this singular point produce a numerical boundary layer about one cell thick along the surface of the step. Woodward and Colella reduce this effect by resetting the cells immediately behind the corner to conserve entropy and the sum of enthalpy and specific kinetic energy through the rarefaction. However, we are less interested here in reproducing the exact solution than in verifying the code and examining the behavior of such numerical effects as resolution is increased, so we do not apply this additional boundary condition. The errors near the corner result in a slight over-expansion of the gas there and a weak oblique shock where this gas flows back toward the step. At all resolutions we also see interactions between the numerical boundary layer and the reflected shocks that appear later in the calculation.

The shock reaches the top wall at $t\approx 0.65$. The point of reflection begins at $x\approx 1.45$ and then moves to the left, reaching $x\approx 0.95$ at $t=1$. As it moves, the angle between the incident shock and the wall increases until $t=1.5$, at which point it exceeds the maximum angle for regular reflection ($40^\circ$ for $\gamma=1.4$) and begins to form a Mach stem. Meanwhile the reflected shock has itself reflected from the top of the step, and here too the point of intersection moves leftward, reaching $x\approx 1.65$ by $t=2$. The second reflection propagates back toward the top of the grid, reaching it at $t=2.5$ and forming a third reflection. By this time in low-resolution runs, we see a second Mach stem forming at the shock reflection from the top of the step; this results from the interaction of the shock with the numerical boundary layer, which causes the angle of incidence to increase faster than in the converged solution. Figure 35.20 compares the density field at $t=4$ as computed by FLASH using several different maximum levels of refinement. Note that the size of the artificial Mach reflection diminishes as resolution improves.

Figure: Density at $t=4$ in the Emery wind tunnel test problem, as computed with FLASH using several different levels of refinement.

Figure: Detail of the Kelvin-Helmholtz instability seen at $t=3$ in the Emery wind tunnel test problem for several different levels of refinement.

The shear cell behind the first (“real”) Mach stem produces another interesting numerical effect, visible at $t \ge 3$ -- Kelvin-Helmholtz amplification of numerical errors generated at the shock intersection. The waves thus generated propagate downstream and are refracted by the second and third reflected shocks. This effect is also seen in the calculations of Woodward and Colella, although their resolution was too low to capture the detailed eddy structure we see. Figure 35.21 shows the detail of this structure at $t=3$ on grids with several different levels of refinement. The effect does not disappear with increasing resolution, for three reasons. First, the instability amplifies numerical errors generated at the shock intersection, no matter how small. Second, PPM captures the slowly moving, nearly vertical Mach stem with only 1-2 cells on any grid, so as it moves from one column of cells to the next, artificial kinks form near the intersection, providing the seed perturbation for the instability. Third, the effect of numerical viscosity, which can diffuse away instabilities on course grids, is greatly reduced at high resolution. This effect can be reduced by using a small amount of extra dissipation to smear out the shock, as discussed by Colella and Woodward (1984). This tendency of physical instabilities to amplify numerical noise vividly demonstrates the need to exercise caution when interpreting features in supposedly converged calculations.

Table 35.6: Runtime parameters used with the WindTunnel test problem.

Variable	Type	Default	Description
sim_pAmbient	real	1	Ambient pressure ($p_0$)
sim_rhoAmbient	real	1.4	Ambient density ($\rho_0$)
sim_windVel	real	3	Inflow velocity ($u_0$)

Finally, we note that in high-resolution runs with FLASH, we also see some Kelvin-Helmholtz roll up at the numerical boundary layer along the top of the step. This is not present in Woodward and Colella's calculation, presumably because their grid resolution was lower (corresponding to two levels of refinement for us) and because of their special treatment of the singular point.

35.1.8 The Shu-Osher problem

The Shu-Osher problem (Shu and Osher, 1989) tests a shock-capturing scheme's ability to resolve small-scale flow features. It gives a good indication of the numerical (artificial) viscosity of a method. Since it is designed to test shock-capturing schemes, the equations of interest are the one-dimensional Euler equations for a single-species perfect gas.

In this problem, a (nominally) Mach 3 shock wave propagates into a sinusoidal density field. As the shock advances, two sets of density features appear behind the shock. One set has the same spatial frequency as the unshocked perturbations, but for the second set, the frequency is doubled. Furthermore, the second set follows more closely behind the shock. None of these features is spurious. The test of the numerical method is to accurately resolve the dynamics and strengths of the oscillations behind the shock.

The shu_osher problem is initialized as follows. On the domain $-4.5 \leq x \leq 4.5$, the shock is at $x=x_{{\tt s}}$ at $t=0.0$. On either side of the shock,

$$\begin{array}{lcccc} & x \leq x_{\tt s} & & & x > x_{\tt s} \... ...& & & p_{\tt R} u(x) & u_{\tt L} & & & u_{\tt R} \end{array}$$ (35.14)

where $a_{\rho}$ is the amplitude and $f_{\rho}$ is the frequency of the density perturbations. The gamma equation of state alternative implementation is used with gamma set to 1.4. The runtime parameters and their default values are listed in Table 35.7. The initial density, $x$-velocity, and pressure distributions are shown in Figure 35.22.

Table 35.7: Runtime parameters used with the shu_osher test problem.

Variable	Type	Default	Description
posn	real	-4.0	Initial shock location ($x_{\tt s}$)
rho_left	real	3.857143	Initial density to the left of the shock ( $\rho_{\rm L}$)
rho_right	real	1.0	Nominal initial density to the right ( $\rho_{\rm R}$)
p_left	real	10.33333	Initial pressure to the left ($p_{\rm L}$)
p_right	real	1.0	Initial pressure to the right ($p_{\rm R}$)
u_left	real	2.629369	Initial velocity to the left ($u_{\rm L}$)
u_right	real	0.0	Initial velocity to the right ($u_{\rm R}$)
a_rho	real	0.2	Amplitude of the density perturbations
f_rho	real	5.0	Frequency of the density perturbations

The problem is strictly one-dimensional; building 2-d or 3-d executables should give the same results along each $x$-direction grid line. For this problem, special boundary conditions are applied. The initial conditions should not change at the boundaries; if they do, errors at the boundaries can contaminate the results. To avoid this possibility, a boundary condition subroutine was written to set the boundary values to their initial values.

Figure: Initial density, $x$-velocity, and pressure for the Shu-Osher problem.

The purpose of the tests is to determine how much resolution, in terms of mesh cells per feature, a particular method requires to accurately represent small scale flow features. Therefore, all computations are carried out on equispaced meshes without adaptive refinement. Solutions are obtained at $t=1.8$. The reference solution, using 3200 mesh cells, is shown in Figure 35.23. This solution was computed using PPM at a CFL number of 0.8. Note the shock located at $x \simeq 2.4$, and the high frequency density oscillations just to the left of the shock. When the grid resolution is insufficient, shock-capturing schemes underpredict the amplitude of these oscillations and may distort their shape.

Figure: Density, $x$-velocity, and pressure for the reference solution at $t=1.8$.

Figure 35.24 shows the density field for the same scheme at 400 mesh cells and at 200 mesh cells. With 400 cells, the amplitudes are only slightly reduced compared to the reference solution; however, the shapes of the oscillations have been distorted. The slopes are steeper and the peaks and troughs are broader, which is the result of overcompression from the contact-steepening part of the PPM algorithm. For the solution on 200 mesh cells, the amplitudes of the high-frequency oscillations are significantly underpredicted.

Figure 35.24: Density fields on 400 and 200 mesh cells from the PPM scheme.

35.1.9 Driven Turbulence StirTurb

The driven turbulence problem StirTurb simulates homogeneous, isotropic and weakly-compressible turbulence. Because theories of turbulence generally assume a steady state, and because turbulence is inherently a dissipative phenomenon, the fluid must be driven to sustain a steady-state. This driving must be done carefully in order to avoid introducing artifacts into the turbulent flow. We use a relatively sophisticated stochastic driving method originally introduced by Eswaran & Pope (1988). The initial conditions sets up a homogeneous background. The resolution used for this test run was $32^3$, and the boundary conditions were periodic. The Table 35.8 shows values the runtime parameters values to control the amount of driving, and the Figures Figure 35.25 and Figure 35.26 show the density and x velocity profile of an xy plane in the center of the domain.

Table 35.8: Runtime parameters used with the Driven Turbulence test problem.

Variable	Type	Value	Description
st_stirmax	real	25.1327	maximum stirring wavenumber
st_stirmin	real	6.2832	minimum stirring wavenumber
st_energy	real	5.E-6	energy input per mode
st_decay	real	0.5	correlation time for driving
st_freq	integer	1	frequency of stirring

Figure 35.25: Density profile for the StirTurb driven turbulence problem.

Figure 35.26: velocity along X dimension for the StirTurb driven turbulence problem.

35.1.10 Relativistic Sod Shock-Tube

The relativistic version of the shock tube problem (RHD_Sod) is a simple one-dimensional setup that involves the decay of an initially discontinuous two fluids into three elementary wave structures: a shock, a contact, and a rarefaction wave. As in Newtonian hydrodynamics case, this type of problem is useful in addressing the ability of the Riemann solver to check the code correctness in evolving such simple elementary waves.

We construct the initial conditions for the relativistic shock tube problem as found in Martí & Müller (2003). We use an ideal equation of state with $\Gamma=5/3$ for this problem and the left and right states are given:

$$\begin{array}{lclcclcl} \rho_{\rm L} &=& 10 & & & \rho_{\rm R... ... L} &=& 40/3 & & & p_{\rm R} &=& 2/3 \times 10^{-6} \end{array}$$ (35.15)

and the computational domain is $0 \le x \le 1$. The initial shock location is at $x=0.5$ (halfway across a box with unit dimensions).

In FLASH, the RHD Sod problem (RHD_Sod) uses the runtime parameters listed in Table 35.9:

Table 35.9: Runtime parameters used with the RHD_Sod test problem.

Variable	Type	Default	Description
sim_rhoLeft	real	10	Initial density to the left of the interface ( $\rho_{\rm L}$)
sim_rhoRight	real	1.0	Initial density to the right ( $\rho_{\rm R}$)
sim_pLeft	real	40/3	Initial pressure to the left ($p_{\rm L}$)
sim_pRight	real	$2/3 \times 10^{-6}$	Initial pressure to the right ($p_{\rm R}$)
sim_uLeft	real	0	Initial velocity (perpendicular to interface) to the left ($u_{\rm L}$)
sim_uRight	real	0	Initial velocity (perpendicular to interface) to the right ($u_{\rm R}$)
sim_xangle	real	0	Angle made by interface normal with the $x$-axis (degrees)
sim_yangle	real	90	Angle made by interface normal with the $y$-axis (degrees)
sim_posn	real	0.5	Point of intersection between the interface plane and the $x$-axis

Figure 35.27, Figure 35.28, and Figure 35.29 show the results of running the RHD Sod problem on a one-dimensional uniform grid of size 400 at simulation time $t=0.36.$ In this run the left-going wave is the rarefaction wave, while two right-going waves are the contact discontinuity and the shock wave. This configuration results in mildly relativistic effects that are mainly thermodynamical in nature.

The differences in the relativistic regime, as compared to Newtonian hydrodynamics, can be seen in a curved velocity profile for the rarefaction wave and the narrow constant state (density shell) in between the shock wave and contact discontinuity. Numerically, it is particularly challenging to resolve the thin narrow density plateau, which is bounded by a leading shock front and a trailing contact discontinuity (see Figure 35.27).

Figure: Density of numerical solution to the relativistic Sod problem at time $t=0.36.$

Figure: Pressure of numerical solution to the relativistic Sod problem at time $t=0.36.$

Figure: Normal velocity of numerical solution to the relativistic Sod problem at time $t=0.36.$

35.1.11 Relativistic Two-dimensional Riemann

The two-dimensional Riemann problem (RHD_Riemann2D), originally proposed by Schulz et al. (1993), involves studying interactions of four basic waves that consist of shocks, rarefactions, and contact discontinuities. The initial condition provided here is based on Migone et al. (2005) producing these elementary waves at every interface. The setup of the problem is given on a rectangular domain of size $[-1,1]\times[-1,1]$, which is divided into four constant state subregions as:

$$(\rho, p, v_x, v_y) = \left\{ \begin{array}{l@{\quad}l} (0.5,1.0,... ...0,0.0) & 0.0 \leq x \leq 1.0,\;\; 0.0 \leq y \leq 1.0 \end{array} \right. ,$$ (35.16)

where $\rho_1=5.477875 \times 10^{-3}$ and $p_1=2.762987 \times 10^{-3}$. An ideal EOS is used with the specific heat ratio $\Gamma=5/3$.

The solution obtained at $t=0.8$ in Figure 35.30 shows that the symmetry of the problem is well maintained. The two shocks are propagated from the upper left and the lower right regions to the upper right region, yielding continuous collisions of shocks at the upper right corner.The curved shock fronts are transmitted and formed in the diagonal region of the domain. The lower left region is bounded by contact discontinuities. By the time $t=0.8$ most of regions are filled with shocked gas, whereas there are still two unperturbed regions in the lower left and upper right regions.

Figure: Log of density of numerical solution to the relativistic 2D Riemann problem at time $t=0.8$. The solution was resolved on AMR grid with 6 levels of refinements

35.1.12 Flow Interactions with Stationary Rigid Body

The stationary rigid body is only implemented and tested in the unsplit hydro solver Sec:unsplit hydro algorithm. It is possible that the unsplit staggered mesh MHD solver Sec:usm_algorithm can support the rigid body but we have not tested yet.

35.1.12.1 NACA Airfoil

Flow simulations over a series of NACA airfoils (or a simple flat plate) can be obtained using the implementation of a stationary rigid body in the unsplit hydro solver described in Sec:StationaryRigidBody. In this example, the cambered NACA2412 airfoil is initialized with positive unity values indicating the part of the domain that belongs to a stationary rigid body. The rest of the domain is assigned negative unity values to indicate it as a flow region for the unsplit hydro solver. At the surface of the rigid body, a reflecting boundary condition is applied in order to represent the fact that there is no flow penetrating the rigid object.

Plots in Figure 35.31(a) - Figure 35.33(b) illusrate Mach number and pressure plots over the airfoil with the three different initial Mach numbers, 0.65, 0.95 and 1.2 at $t=1.8.$ By this time, the flow conditions have reached their steady states. For Mach number = 0.65, the critical Mach number has not yet been obtained and the flow over the airfoil is all subsonic as shown in Figure 35.31(a). Since the airfoil is asymmetric and cambered, we see there are pressure gradients across the top and bottom surfaces even at zero angle of attack in Figure 35.31(b). These pressure gradients (higher pressure at the bottom than the top) generate a lift force by which an airplane can fly defying gravity.

At Mach number reaching 0.95 as shown in Figure 35.32(a) there are local points that are supersonic. This indicates that the critical Mach number for the airfoil is between 0.65 and 0.95. In fact, one can show that the critical Mach number is around 0.7 for the NACA2412 airfoil. We see that there is a development of a bow shock formation in front of the airfoil. A formation of a subsonic region between the bow shock and the nose of the airfoil is visible in the Mach number plot. Inside the bow shock, a sonic line at which the local flow speed becomes the sound speed makes an oval shape together with the bow shock. In both the Mach number and pressure plots, a strong wake forms starting from the top and bottom of the surfaces near the trailing edge. The wake is hardly visible for Mach number 0.65 in Figure 35.31(a) and Figure 35.31(b). Normal shock waves have formed steming from the trailing edge as seen in Figure 35.32(a) and Figure 35.32(b).

At Mach number 1.2 the flow becomes supersonic everywhere. In this case, the shape of the bow shock becomes narrower and there are much larger supersonic pockets developed on the top and bottom surfaces with a smaller subsonic region between the bow shock and the nose region.

Figure: NACA2412 in Mach number 0.65 flow at 0 degree angle of attack problem at $t=1.8$ (a) Mach number (b) Pressure

[a] [b]

Figure: NACA2412 in Mach number 0.95 flow at 0 degree angle of attack problem at $t=1.8$ (a) Mach number (b) Pressure

[a] [b]

Figure: NACA2412 in Mach number 1.2 flow at 0 degree angle of attack problem at $t=1.8$ (a) Mach number (b) Pressure

[a] [b]

35.1.12.2 Solid Objects in Sedov Explosion

Another problem for testing a stationary rigid body in a simulation is to consider the Sedov-like explosion in a chamber surrounded by a solid wall with holes The wall is shown as red blocks with white boundry in Figure 35.34(a) and Figure 35.34(b). The simulation was done on a uniform Cartesian grid with 300 cells on each direction. Three holes in the wall subdivide it into four different stationary solid bodies in a square computational domain $[0,1.5] \times [0,1.5]$. The explosion goes off at the origin and generate shock waves inside the chamber. In later time, when the shock waves pass though the three holes in the wall, turbulence effects are triggered from the interaction between the fluid and the wall and enhance vortical fluid motions.

One important thing in this problem is to keep the given symmetry throughout the simulation. The flow symmetry across the diagonal direction is well preseved in Figure 35.34(a) and Figure 35.34(b).

Figure: The Sedov explosion in a chamber surrounded by a wall with holes. (a) Density plot at $t=0.0$ sec. (b) Denstiy plot at $t=0.1$ sec.

[a] [b]