To test the striping instability in codes with exact geometry, requires a test where the shock traversing the grid is nearly perfectly aligned, and infact is curved with respect to the coordinate grid and hence is intrinsicly 2-D in nature.
To this end a new striping test is proposed here, using an adiabatic wall shock as the test model. A low density fast shock is reflected to the left, off a dense layer on the right. The inflow from the left that generates the shock is steady, but with some variations that result in a slightly non-planar shock.
The inflow is modified so that in the cetral rows and at the top and bottom of the grid in a periodic fashion, the inflow is between 5 and 10 percent faster over three cell rows. This produces 'dimples' in the resulting reflected shock and the shock forn itself is slightly non-planar. The enhanced velocity rows are described with a parameter, delat_v, the fractional increase in row speed over the main grid inflow speed. Only small velocity variations are needed to induce strong striping.
(Note this setup is similar to the radiative test in Sutherland et al 2003, but is adiabatic, dimensionless and hence hopefully more reproducable by other codes.)
- Adiabatic Index. Gamma = 5/3
- Domain: 0.0 < x <2.0, 0.0 < y < 1.0
- 288 x 144 cells
- Density: 1.0
- Pressure: 1.0
- Velocity: v_x = v_in = 15.0, v_y = 0.0, modified below
- Y boundaries are periodic, X boundaries are fixed, constant values
- Density: 300.0
- Pressure: 300.0
- Velocity: v_x = 0.0, v_y = 0.0
Initial region boundary
- x_0 = 1.8
Inflow velocity, left, x= 0.0 boundary
- v_in = 15.0
The velocity on the left, x = 0.0, boundary, and throughout the Left Region initially, v_x = v_in, is modified on specific grid rows as follows, taking j = 0... (n_y-1), where n_y is, the number of y cells in grid, (here n_y = 144 for the standard test. The parameter delta_v = 0.1 for the standard test, and is the fractional increase in velocity in some rows on the grid. v_y = 0.0 everywhere.
- v_x[j] = v_in * (1.0 + delta_v), for (j mod (n_y/2)) == 0, ie j = 0, 72
- v_x[j] = v_in * (1.0 + 0.5*delta_v), for ((j-1) mod (n_y/2)) == 0, ie j = 1, 73
- v_x[j] = v_in * (1.0 + 0.5*delta_v), for ((j+1) mod (n_y/2)) == 0, ie j = 71, 143
- v_x[j] = v_in everywhere else
Right, x = 2.0, boundary
- v_x = 0.0, v_y = 0.0 everywhere
- CFL number: 0.8, (Standard)
- Flattening: minimum 0.0, maximum 1.0 (standard)
- Alpha: 0.0, 1e-4, 0.001, 0.01
- Time, t = 0.3
The simulation is run with all striping control code turned off, then on, using a smoothing constant of Alpha = 0.0001, 0.001 and 0.01.
With the numerical diffusion used in Fyris Alpha, the parameter alpha is the fraction of a cell quantity value that is split off, divided in two, and added to the two neighbouring cells, conservatively smoothing the quantity. The smoothing is only applied in regions where the divergence of the orthogonal velocity field is great enough, owing to a shock orthoganal to the direction of the code sweep (see Sutherland et al 2003). The improvement here over the 2003 algorithm is to use alpha times the magnitude of the velocity divergence for smoothing, rather than a fixed alpha, above a given threshold, to simulate a viscosity. Additionally Fyris Alpha only applies the smoothing to the density and velocity variables, and does not to diffuse the pressure. In very strong shocks it is errors in the large post shock pressures that may be driving the striping instability, so the propagation of large pressure differences is limited here by not diffusing the pressure variable.
In addition, the velocity of the contact discontinuities, derived from the cell boudnary Riemann problems, is set to 0.0 in sweeps of cells orthogonal to shock fronts, again simulating high transverse viscosity in the shock front regions and reducing the buildup of stripes. In practice this is similar to the H-Correction ofSanders et al. (1998), but applied in a split sweep semi-lagrangian context. The H-Correction in addition to the alpha-viscosity dissipation in shock fronts, allows alpha to be much smaller, 0.001 to 0.01 in Fyris, compared to 0.075 used routinely in ppmlr.
Note, that as the smoothing and viscosity is confined to orthogonal sweeps of strong shock fronts, the majority of the grid is free of any additional disspation, and the outcome, once striping is suppressed, is more or less independent of the amount of dissipation introduced. The results for alpha = 0.001 and 0.01 are more or less identical with only minor levels of extra diffusion visible for alpha = 0.01.
Because this smoothing operator is applied per cycle of computation, it is not totally independent of the courant number or the absolute strength of the shocks being stabilised. In the case of the 2-D Noh shock, with a Mach number of around 1000, a value of alpha = 0.05 was needed to prevent the circular shock from breaking up under the striping instability. The striping is very strong for this essentially infinite strenght shock problem, and it almost certainly the origin of the failure for many of the L&W test codes to compute the 2-D Noh problem correctly.
No striping control
Striping control, alpha = 1.0e-4
Striping control, alpha = 0.001
Striping control, alpha = 0.01