Linear Waves

References

Plasma Physics by P.A. Sturrock, Chapter 14, Cambridge U Press, 1994 The Physics of Astrophysics, Volume II: Gas Dynamics by F. Shu, Chapter 22, University Science Books, 1992

Description

This test problem is the propagation of linear waves in a background medium. It is initiated by choosing from any of the wave modes for ideal hydro/MHD (e.g., sound waves, Alfvén waves, magnetosonic waves) and perturbing a uniform medium to produce the desired wave. While these waves are only linear in amplitude, preventing us from testing the algorithm in the nonlinear regime, this test provides a good measure of the convergence rate of the code. One can compare the numerical solution after some integer number of wave periods with the initial condition. Differences from the initial condition will show up from dispersion error and/or diffusion error. Both of these errors should converge at the order of the algorithm (for a second-order algorithm, we should have a convergence rate of -2). One can also compare the L1 norm and convergence rates for two waves of the same type but propagating in opposite directions. Since numerical error is symmetric in Athena3D, the L1 norm should be identical between two such waves. Correct convergence is a necessary (though, not sufficient) condition with which to test an algorithm.

Set up

Domain 1D: 0.0 ≤ x ≤ 2.236068 2D: 0.0 ≤ x ≤ 2.236068, 0.0 ≤ y ≤ 1.118034 3D: 0.0 ≤ x ≤ 2.236068, 0.0 ≤ y ≤ 1.118034, 0.0 ≤ z ≤ 0.745356 Boundary conditions Periodic everywhere Equation of state The linear waves problem can be run with an isothermal or adiabatic equation of state. For adiabatic runs, we use γ = 5/3. Initial density ρ = 1 everywhere with a perturbation added (see "Perturbations" section below). Initial pressure For an adiabatic gas, P = 1/γ with a perturbation added (see below). These values for P and ρ give us an adiabatic sound speed of 1. For an isothermal equation of state, the value of P doesn't need to be set in Athena3D. In this case, we choose the isothermal sound speed to be 1. Initial velocity To produce a pure wave mode without any background flow, we set all velocities to zero, with appropriate perturbations added (see below). One can set the background velocity to the negative of the wave speed to produce a standing wave. For the contact waves (where there is only advection), the background velocity is set to 1. MHD Components This test problem can be done with or without MHD components. Bx = 1 with perturbation added (see below). By = 21/2 with perturbation added (see below). Bz = 0.5 with perturbation added (see below). Perturbations To initialize the proper wave mode, a sinusoidal perturbation is added to the background state described by the above initial conditions. We describe this perturbation in 1D here. The extension to 2D and 3D is presented in the "Rotated Waves" section. The perturbation can be written in array form as, δU = A R sin[(2π/λ)x], with A being the perturbation amplitude, which we choose to be 1x10-6. The wavelength of the perturbation is given by λ and in 1D, λ = Lx, where Lx is the length of the domain in the x direction. The array δU is a perturbation to the array of conserved variables, U = [ρ, ρ vx, ρvy, ρvz, E, Bx, By, Bz]. The appropriate components of this array are dropped depending on what physics the code is run with (e.g., for isothermal gas, the total energy, E, component is dropped). R is the right eigenvector corresponding to the wave mode chosen. In other words, R contains the relative perturbation amplitudes and signs for each of the variables in the conserved variable array. In Athena3D, we can obtain R directly from inputting the background state values into the linearized Riemann solver, which then gives us the eigenvalues and eigenvectors for the desired wave modes. Rotated Waves We run the linear waves problem in both 2D and 3D in addition to 1D. For the 2D problem, we take the 1D configuration described above and rotate it by an angle α with respect to the x-axis, where α = tan-1(Lx/Ly). The wavelength of the perturbation is then λ = Lxcos(α). These parameters were chosen so that the rotated wave is periodic in x and y. The rotation by α is applied to the components of the momentum and magnetic field. To guarantee that the magnetic field has zero divergence initially, we calculate the z component to a vector potential from the rotated magnetic field components. This potential is located at the zone corners so that the finite difference form of the curl operation determines the field components at the cell interfaces. These components are the fundamental magnetic field variables in Athena3D. For the 3D problem, we apply the same rotation by α with respect to the x-axis and an additional rotation by β = tan-1[(Lxcos(α)+Lysin(α))/(2Lz)] with respect to the xy plane. The wavelength of the perturbation is λ = Lxcos(α)cos(β). We calculate all of the components of a vector potential (more than just the z component is needed in 3D) from the rotated magnetic field components, and then use this vector potential to guarantee the magnetic field has zero initial divergence. Convergence Testing As discussed above, the primary reason for performing this test is to measure the convergence of the algorithm. In our simulations, we evolved the waves for 4 wave periods, after which we calculate the quantity ΔUijk = |Uijk(t = 4 periods) - Uijk(t = 0)|, which is the absolute value of the difference between the array of conserved variables after 4 wave periods and the initial conserved variables. The ijk subscripts are indices denoting the i-th grid cell in the x direction, j-th grid cell in the y direction, and the k-th grid cell in the z direction. We then average this difference over the entire grid and add the components of the conserved variable array in quadrature, giving us our measure of the L1 error: L1 norm = [∑m(∑ijkΔUijk,m)2]1/2/(NxNyNz), where m indicates the m-th component of the ΔUijk array, and Nx, Ny, and Nz are the number of grid points in the x, y, and z directions, respectively. Note that all waves shown below are left moving wave modes, except for the contact modes.

Results

1D Isothermal Hydro We ran Athena3D with an isothermal equation of state and no magnetic field evolution in 1D at resolutions Nx = 100, 200, 400, 800, and 1600. Plotted below is the L1 norm versus x resolution (Nx) for a sound wave and a vy contact wave (solid lines with squares). The dashed line shows the slope of ideal -2 convergence. The convergence is close to or better than ideal, with the exception of the highest resolution simulations of the sound wave. For these resolutions, the error appears to level off with increasing resolution. This effect is due to roundoff error beginning to dominate over truncation error (see "Summary" for more details).

1D Adiabatic Hydro Plotted below is the convergence rate for the same wave modes and resolutions as above, but with an adiabatic equation of state. Again, the dashed line shows the slope of -2 convergence. The convergence rates appear very similar to the isothermal case. The sound wave has slightly larger errors in the adiabatic case than in the isothermal case. The vy contact wave is identical between the two cases, which makes sense since the equation of state shouldn't affect contact modes.

1D Isothermal MHD Here are the convergence rates for the 1D fast magnetosonic, Alfvén, and slow magentosonic waves in an isothermal gas for the same resolutions as the previous simulations.

1D Adiabatic MHD Here are the convergence rates for the same waves as the 1D isothermal MHD simulations, but with an adiabatic equation of state and with a density contact wave added.

To show the effect of numerical resolution on the actual wave itself, we have plotted the amplitude of the x velocity versus x for a fast magnetosonic wave in 1D. The plot shows the wave after four wave periods for several resolutions. The resolutions used for the convergence plots above are actually too large (i.e., the wave is really well resolved) to clearly show the effect of changing resolution. Therefore, we run the fast wave for resolutions Nx = 8 (black, dashed line), 16 (black line), 32 (red line), 64 (green line), and 128 (blue line). The plot shows the primary source of numerical error is diffusion error. By Nx = 64, the wave is pretty well resolved.

2D Adiabatic Hydro We have run a sound wave and a density contact wave in 2D hydro with an adiabatic equation of state. The lowest resolution simulation is Nx = 32, Ny = 16, and we increment the resolution by a factor of two in each dimension up to Nx = 512, Ny = 256 to get the convergence plots below. The geometry and resolutions guarantee that the wave propagates along the diagonal of the box and that the fluxes returned by the Riemann solver are different in the x and y directions.

2D Adiabatic MHD Here are the same MHD waves as in the 1D case, but run in 2D at the same resolutions used in the 2D adiabatic hydro simulations.

3D Adiabatic Hydro Here again is the adiabatic hydrodynamic sound and density contact waves, but this time run in 3D. The lowest resolution is Nx = 24, Ny = 12, Nz = 8, and we increment the resolution by a factor of two in each dimension up to Nx = 192, Ny = 96, Nz = 64. Again, the geometry and resolutions guarantee that the fluxes in the x, y, and z directions will all be different.

3D Adiabatic MHD Here are the same MHD waves as in the 1D case, but run in 3D at the same resolutions used in the 3D adiabatic hydro simulations.

Summary All linear wave simulations show close to -2 convergence. The 1D waves show very rapid convergence up to the highest resolutions, where roundoff error begins to dominate over truncation error. The dominance of roundoff error gives an indication of how small the truncation error can be in Athena3D (these simulations were run in double precision). The dominate source of truncation error appears to be diffusion (at least in 1D).

Return to Test Suite page.