Test simulations : MHD/HD with Godunov, TVD-MUSCL etc.

• Orszag-Tang vortex test problem [Orszag and Tang, J. Fluid Mech. (1979)], 2D-MHD, with TVD-MUSCL method

  mpeg file (1.1MB)

  Snapshot at time ∼ π : The animation movie shows temperature for time 0 to 10 simulated with 256x256 grids covering 2π-by-2π space. Mind that color is controlled by the maximum and minimum in each frame, thus does not represent absolute values.

  The small-scale discontinuities and vortices, especially at chaotic regime t > 5, are really sensitive to choices of mesh size and computation scheme. The images below show the temperature at t = 7 simulated with, from the left, 64x64, 128x128, 256x256, 512x512 grids, respectively. The identical code, TVD-MUSCL with Powell correction of divergence B, was used.
  
  
• MHD aligned rotor [c.f. Balsara and Spicer, JCP vol.149 pp.270-292 / vol.153 p.671 (1999)], 2D-MHD, with TVD-MUSCL method.

  Left plot and the background image of this page show the snapshot of temperature at t ∼ 1 and 1.5, respectively.

  
  
• Interaction between oblique MHD shock and dense cloud [c.f. Dai and Woodward, ApJ vol.494 p.317 (1998)], 2D-MHD, with TVD-MUSCL method

  Temperature at t ~ 0.055.

  
  
• Kelvin-Helmholtz Instability, 2D-MHD, with TVD-MUSCL method
  The hydrodynamics case is very chaotic, as seen in this movie (mpeg ; 1.7MB). The below are snapshots of temperature at t = 0, 1/2, 1, 2, 5. Here a small mass density enhancement, instead of small random flow, is given at the lower border to induce/seed instability.
  Under the presence of the magnetic field, results are quite different. Plots, show the simulated temperature at t = 0, 1 and 4 (with 256x256 grids) of the weakly-magnetized plasma. A uniform horizontal field (beta-ratio is 20) is given as the initial field. The initial plasma (t = 0) is identical to the HD case above. This movie (mpeg ; 655kB) shows temperature at t = 0 ∼ 10. Magnetic field works to stabilize system, suppressing the evolution of small-scale instabilities.
  
  
• Rayleigh-Taylor instability, 2D-HD/MHD, with TVD MUSCL method

  HD case mpeg file (679kB)

  Vertical mpeg file (191kB)

  Horizontal mpeg file (117kB)

  Oblique mpeg file (151kB)

  From left shown are the snapshots at t ∼ 12.75 of hydrodynamics case (without magnetic field), and three magnetohydrodynamics cases starting with vertical, horizontal and oblique magnetic field. The initial magnetic field in all MHD cases is uniform and, its strength is set so that the beta ratio will be 250 at the boundary between the high- and low-density regimes. All four cases start with the identical initial plasma distribution with tiny "single"-mode symmetric fluctuation of bulk flow. The movies in the left show the temporal evolutions from t = 0 to 15. Here again, (tension of) magnetic field works to suppress hydrodynamical instability with different ways depending on the direction of field.

  
  

  mpeg file (481kB)

  When the initial magnetic field has only the component perpendicular to the plotted plane, the magnetic field will behave as if it is a plasma with specific heat ratio of 2. In this case, no magnetic tension will arise, and thus instability will not be suppressed.

  
  

  HD case mpeg file (995kB)

  MHD case mpeg file (429kB)

  The evolutions of RT instability of HD (left) and MHD (right) cases starting with multiple-wavelength small-amplitude fluctuation of plasma : The MHD case starts with uniform horizontal magnetic field (beta-ratio is about 1000). Note that the difference in file size of mpeg files well represents the differences in degree of complexity in these two cases.

  
  
• 1-D MHD shock/discontinuity problem with TVD-MUSCL method [c.f. Ryu and Jones, ApJ vol.442 p.228-258 (1995)],

  mpeg file (1.4MB)

  The initial values at left and right side, in a order [Rho, Vx, Vy, Vz, Bx, By, Bz, Pg], are [1, 0, 0, 0, 1.3, 1.0, 0, 1.0] and [0.4, 0, 0, 0, 1.3, -1, 0, 0.4], respectively. The specific heat ratio (γ) is set to be 5/3. The non-reflection boundary treatment based on the method of characteristics is applied at both ends.

Other test simulations

• Non-reflection boundary treatment by means of the projected normal characteristics method [c.f. Nakagawa, (1980); Wu et al., (1983); Han et al., (1988)]

  mpeg file (309kB)

  mpeg file (264kB)

  The two plots show the snapshot of the gas pressure; the MHD (left) and HD(right) waves propagate, responding to the impulse given at the center. The non-reflection boundary treatment based on the "projected normal characteristics method" is implemented on the four boundaries. Mainly because this method uses only the normal component of the characteristic curves for simplicity, a perfect non-reflect is not obtained. However, most wave (energy) goes outside the computation domain with minor reflection, and such small reflection waves vanish after reflecting twice or so.

  In this way, this method can provide MHD/HD solution simultaneously satisfying the governing (hyperbolic) equations and the concept of the characteristics on the "physical boundary". This method is very helpful for the simulation studies of the solar corona in which we usually assume the photosphere (or "coronal base" at some height) as one of the computational boundaries.
  
  
• WENO (weighted essentially non oscillatory) method, of 5th-order [c.f. Jiang and Wu, JCP vol. 150-2 pp.561-594(1999)]

  The snapshots of Orszag-Tang vortex problem at time ∼ π.    The left plot is a snapshot of the simulated temperature with the 5th-order weighted essentially non-oscillatory (WENO) method and the Roe's method. The right one is obtained with 3rd-order TVD-MUSCL (here, with "minmod" function) and thus identical to the plot at the top part of this page.

     Both schemes work pretty well for solving many kinds of MHD/HD problems, and there seems less differences between the two plots. Looking at details, though, we can find some.
     With the TVD method, the nonlinear steepening at shock discontinuity is better treated, and the artificial viscosity to stabilize the computation can be minimized. The WENO and its earlier version, ENO, can better treat the maxima and minima of smoothed curves like sinusoidal without generating false extrema. The WENO-Roe is a combined version of these two methods, picking up virtues of each.
  
  
• HLL (Harten-Lax-van Leer) approximate Riemann solver.
  ... and HLLC (HLL contact), HLLCD (HLL discontinuities)
  coming soon ...
• CT (constraint trasnport) method
  ... for solving divergence-free field governed by induction equation
  coming soon ...
• Conjugate Gradient method(s)
  coming soon ...

Other simulations : … just for fun

• Particle motions of the two-dimensional hydrodynamical flow with an asymmetric heat source, with two-step Lax-Wendroff method and Lapidus-type artificial viscosity

  mpeg file (978kB)

  An asymmetric heat source is given on the bottom boundary to induce asymmetric motions.    Lax-Wendroff method is numerically diffusive, though, has been long used because it is much less expensive. Computations with this method is generally 4 ∼ 6 times faster than that with the TVD method of Godunov type approximate Riemann solver. This is an advantage for various simulation problems.

  
  
• Supersonic flow (M=10) injection, 2D-HD, with Rusanov/Lax-Friedrichs method

  mpeg file (395kB)

  A tiny density enhancement is given near the injection site so that symmetry will be broken. Lax-Friedrichs method is a simplified version of TVD method. This is a quick and enough accurate scheme.

  
  
• Collision of two fluid bodies, 2D-HD, Godunov TVD-MUSCL

  mpeg file (239kB)

  Just for fun.    Collision of two objects is simulated with compressive fluid code. With the initial pressure being uniform everywhere in the simulation box, the objects are featured by different density (thus, different temperature) values. The ratio of initial densities at the core (inner part) of the left object, the outer part of the left object, the right object, and the surrounding space, is 20 : 10 : 15 : 1. The right object is moving leftward to collide onto the left one. Periodic boundary is applied on all four boundaries. No gravity is taken into account. This simple setting looks working, anyway.