- Email: [email protected]
Application of a total variation diminishing scheme to electromagnetic hybrid particle-in-cell plasma simulation
Computer Physics Communications 183 (2012) 2027–2034
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc
Application of a total variation diminishing scheme to electromagnetic hybrid particle-in-cell plasma simulation Masaharu Matsumoto a,∗ , Yoshihiro Kajimura b , Hideyuki Usui a , Ikkoh Funaki b , Iku Shinohara b a
Department of Computational Science, Kobe University, Japan
b
Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, Sagamihara, Japan
article
info
Article history: Received 18 October 2011 Received in revised form 13 April 2012 Accepted 24 April 2012 Available online 3 May 2012 Keywords: Total variation diminishing Particle-in-cell code Electromagnetic hybrid plasma code Space plasma physics
abstract A discretization procedure for a total variation diminishing (TVD) scheme is introduced to an electromagnetic hybrid particle-in-cell (PIC) plasma simulation code in order to improve the numerical stability and resolution when calculating the plasma flow field in which magnetic field discontinuities (for example, Rankine–Hugoniot jump conditions for shock waves) are generated. In the hybrid PIC code used in this study, ions are treated as particles and electrons are assumed to be an inertia-less (massless) fluid. In the numerical results of one-dimensional test simulations, the TVD scheme significantly prevents non-physical, numerical oscillations, which would ordinarily be produced in the solution when the convection term of the magnetic induction equation in the hybrid PIC code is discretized by central difference schemes at magnetic field discontinuities. Furthermore, a two-dimensional simulation of the global structure of a collision-less bow shock, which is suitable for practical use, makes it possible to clearly capture the bow shock by using the hybrid PIC code with the TVD scheme. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Hybrid particle-in-cell (PIC) plasma simulation models [1,2], in which the various components of a plasma are treated in different ways and the plasma can be coupled to the electromagnetic field, have been widely utilized in plasma science and engineering over the past three decades. In the most common type of hybrid PIC code, the ions are treated kinetically as particles, the electrons are modeled as an inertia-less (mass-less) fluid, and the electromagnetic field is applied to the low-frequency (Darwin) approximation. Such models are designed to reveal the plasma phenomena that occur on shorter time and space scales than can be described by magnetohydrodynamics (MHD) models, although the detailed physical processes on electron inertial scales (for example, Debye length and electron Rarmar radius scales, and electron plasma frequency and electron cyclotron frequency time scales) are not resolved. In other words, hybrid PIC codes focus mainly on the plasma behaviors that occur on ion inertial scales (for example, ion Rarmar radius scales and ion cyclotron frequency time scales), namely, between the MHD and electron inertial scales. These models are used to reasonably simulate plasma phenomena
∗ Correspondence to: Department of Computational Science, Graduate School of System Infomatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan. Tel.: +81 0 78 803 6210; fax: +81 0 78 803 6210. E-mail addresses: [email protected] (M. Matsumoto), [email protected] (H. Usui). 0010-4655/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2012.04.021
in which ion kinetic effects are important. In hybrid PIC codes, the momentum equations for ion particles, generalized Ohm’s law (or momentum equation for an electron fluid), energy equation for an electron fluid, and Faraday’s induction equation (or induction equation for a magnetic field) may be solved numerically. Numerical algorithms and implementations of hybrid PIC codes have been developed and may be divided into several methods: the predictor–corrector scheme suggested by Harned [3], the model suggested by Terasawa et al. [4], the moment method suggested by Winske et al. [5], the QN3D method suggested by Horowitz et al. [6], the current advance method and cyclic leapfrog (CAM-CL) suggested by Matthews [7], amongst others. The predictor–corrector scheme [3] continues to be commonly utilized, even though it was first proposed as a hybrid PIC code. This is because it provides good energy conservation and robustness. However, it is well known that a considerable amount of shortwavelength Whistler noise is generated in this scheme. In order to avoid such noise, subcycling time steps are introduced and different time increments are used for the particles and fields in Terasawa’s model [4]. In other words, the time increment used for advancing the field, which is smaller than that for advancing the particles, is set to a value that is small enough to satisfy the requirement of the characteristic velocity of the Whistler wave, and the field calculations iterate until the particles are advanced. A second-order rational Runge–Kutta method is used to advance both the fields and the particles in this method. The moment method [5] has been proposed as a variation of Terasawa’s model.
2028
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
There are also subcycling time steps, but with a fourth-order Runge–Kutta method. Here, the ion fluid velocity is advanced using an MHD momentum equation. Instead of using the MHD momentum equation for advancing the ion fluid velocity, a fourthorder Adams–Bashforth extrapolation is also used [1]. In the QN3D method [6], functional iteration for advancing the fields is introduced instead of the subcycling time steps. In the CAM-CL scheme [7] that has becomes popular in recent years, the ion fluid velocity is calculated by advancing an extra half time step using a mixed level evaluation of the electric field. Additionally, with the rapid advancement of computer power, hybrid PIC codes have recently been developed that are suitable for massively parallel computation [8] or adaptive mesh refinement [9]. Although many variations of hybrid PIC codes obviously exist, with the exception of the methods mentioned above, the most notable differences between these codes involve the method used to advance the time steps or the methods used for the time integration of the particles and fields. In contrast, insufficient attention has been paid to the discretization of the spatial derivatives that appear in the magnetic induction equation, compared to that of the time derivative. The magnetic induction equation is a type of nonlinear hyperbolic partial differential equation, which has the convective derivative term shown below.
∂ b ∂ f (b) + = 0, ∂t ∂x
(1)
where b and f indicate an arbitrary physical quantity and function, respectively. Unfortunately, as is well known, the numerical analysis of such a hyperbolic conservation equation is not amenable to straightforward time and space integration and may give rise to numerical instabilities and non-physical oscillations if it is integrated in a naive way. In particular, almost all of the standard spatial discretization methods yield numerical results with excessive non-physical oscillations when a discontinuous surface appears in the physical quantity, for example, the Rankine–Hugoniot jump conditions of a shock wave. Therefore, this type of equation requires special care to be integrated numerically. One way to overcome this problem is to introduce an artificial viscosity as a smoothing filter. The idea that stable computations can be accomplished through an artificial viscosity or dissipater was used in an early scheme of computational fluid dynamics (CFD). The artificial viscosity smoothes out the strong gradients and attenuates the spurious oscillation in the solution. However, the difficulty with the artificial viscosity approach is that it is hard to determine an appropriate form that introduces just enough dissipation without causing unnecessary smearing. Over the past two decades, in CFD (including MHD), total variation diminishing (TVD) schemes [10,11], which were first introduced by Harten [12], have been widely used. The TVD schemes seldom generate spurious oscillation for nonlinear hyperbolic conservation laws because these guarantee that the total variation of the numerical solution will not increase under the numerical process in time. Any numerical scheme that gives rise to such oscillations does not satisfy the TVD conditions. For example, the central difference schemes that have been widely used in the convection terms of induction equations in conventional hybrid PIC codes are not TVD schemes. Based on the above, the purpose of this paper is to improve the numerical stability of a hybrid PIC code by applying a TVD scheme to the induction equation. 2. Numerical implementations The intended plasma in this paper is a singly ionized, fully ionized, and collision-less plasma. Fig. 1 shows a timechart of the hybrid PIC code used in this study, which is almost the same as Terasawa’s scheme [4]. Here, X# , V# , N , V, E, and B
Fig. 1. Timechart for the algorithm of the hybrid PIC code used in this paper.
indicate the position of the ion particles, velocity of the particles, number density, flow velocity, electric field, and magnetic field, respectively. Subscript i indicates the physical value of the ions. The general layout of the algorithm for the hybrid PIC code used in this study is as follows. Because the details of this hybrid PIC code can be found in Refs. [1,2], only a brief description is given below. We begin with the positions of the ion particles at time n, their velocities at time n + 1/2, and the electric and magnetic fields at time n. In step 1, the positions of the ion particles are updated from time n to n + 1/2 and then to n + 1. In step 2, the ion density and fluid velocity are calculated from the positions and velocities of the ion particles using the PIC method. In step 3, the magnetic field is updated by solving a magnetic induction equation while the ion density and ion fluid velocity are fixed in time, and then the electric field is calculated from the updated magnetic field. Here, subcycling time steps are used to update the magnetic field. Finally, in step 4, the velocities of the ion particles are updated using the Boris algorithm. Note that, in this paper, we focus exclusively on solving the induction equation in the hybrid PIC code. Therefore, the algorithm of step 3 will mainly be discussed later. For details about the numerical algorithms used in the hybrid PIC code, please refer to Refs. [1,2]. If the Darwin approximation and quasi-neutrality can be adopted for the intended plasma, the momentum equation for a mass-less electron fluid (generalized Ohm’s law), Maxwell’s equations, and the charge neutrality equation are shown as follows. E + Ve × B +
∇ Pe 2Ne
∂B = −∇ × E ∂t j=∇ ×B ∇ ·B=0 N = Ni = Ne ,
=0
(2) (3) (4) (5) (6)
where N , B, V, E, j, and Pe are the number density, magnetic field, flow velocity, electric field, current density, and electron pressure, respectively. Subscripts e represent electrons. Furthermore, the following normalizations are used: the number density, N, is normalized by the initial uniform number density, Nc ; the magnetic field, B, is normalized by the constant field magnitude, Bc ; the velocity, V, is normalized by the characteristic Alfven velocity, VA (=Bc /[µ0 mi Nc ]1/2 , where µ0 and mi are the vacuum permeability and ion mass, respectively); the electric field, E, is normalized by Bc VA ; the current, j, is normalized by eNc VA (e is the electric charge); and the pressure, P, is normalized by the magnetic pressure, B2c /2µ0 . The units for time and space are the reciprocal
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
of the ion cyclotron frequency, ωci−1 , and the ion inertial length, VA /ωci , respectively. The magnetic induction equation is given by substituting Eq. (2) into Eq. (3):
∂B ∇ Pe 1 = ∇ × Vi × B + − (∇ × B) × B . ∂t 2N N
∂Q ∂F + =0 ∂ t ∂ x By Vx By − Vy Bx Q= , F= . Bz Vx Bz − Vz Bx
(7)
1t n (F˜ − F˜jn−1/2 ). 1x j+1/2
(8)
(9)
Let Qjn be the numerical approximation of By or Bz at x = j1x and t = n1t, where 1x is the spatial mesh size and 1t is the time step. F˜ is a numerical flux vector function. The numerical flux for the second-order upwind TVD scheme suggested by Harten and Yee can be written as F˜j+1/2 =
1 2
(Fj+1 + Fj + φj+1/2 )
(10)
φj+1/2 = σ (cj+1/2 ) (gj + gj+1 ) − ψ(cj+1/2 + γj+1/2 )∆j+1/2
(11)
gj = minmod (∆j+1/2 , ∆j−1/2 )
(12)
minmod (x, y) = sgn(x) · max{0, min[|x|, y sgn(x)]}
(13)
γj+1/2
gj+1 − gj = σ (cj+1/2 ) ∆ 0 j+1/2
σ (z ) = ψ(z ) =
1 2
1t 2 ψ(z ) − z 1x
|z | (z 2 + δ 2 )/2δ
∆j+1/2 ̸= 0
(14)
∆j+1/2 = 0
|z | ≥ δ |z | < δ.
Quantity
Value
Time step Grid size Particles Grid number Ion density Magnetic field
0.01 0.5 36/cell 1024 1.0 1.0
schemes. Because the details of TVD schemes can be found in Ref. [10], only a brief description is given in this paper. Extensions to multi-dimensional equations can be obtained using a method similar to that mentioned above. The remaining second and third terms on the right-hand side of the magnetic induction equation (Eq. (7)) are discretized by a second-order central difference scheme. Subcycling time steps are also used [1,2,4]. That is, different time increments are used for update of the particles and fields. 3. Results
Here, note that the equation for Bx is ignored. We will discretize Eq. (8) while the ion fluid velocities (Vx , Vy , Vz ) are fixed in time. Generally, the numerical solvers for a nonlinear system equation of the hyperbolic conservation law are complicated, because there are characteristic velocities based on the number of equations. However, the system equation of Eq. (8) can be treated as two independent scalar equations for By and Bz , because the characteristic velocities of Eq. (8) are multiple roots of Vx . Many existing conservative schemes for the system of Eq. (8) use forward Euler time discretization. For simplicity, this time discretization is used for the discussion in this paper, and the scheme can be written as Qjn+1 = Qjn −
Table 1 Non-dimensional parameters for wave dispersion relation analysis.
Here, Eqs. (4) and (6) and the definition for the current, j = Ni Vi − Ne Ve , are substituted into Eq. (2). The first term on the right-hand side in Eq. (7) represents the convection term, while the second and third terms represent the electron pressure gradient term and Hall term, respectively. Note that the sum of the left-hand side and the first term on the right-hand side is the same hyperbolic form as Eq. (1). We attempt to apply the TVD scheme to the convection term of the magnetic induction equation. For simplicity, a onedimensional induction equation in which the time derivative and convection terms are extracted is considered. Then, the nonlinear system equation for the magnetic induction equations can be obtained, as shown below.
2029
(15)
(16)
Here, cj+1/2 indicates the characteristic velocity. In this case, cj+1/2 is equal to Vxj+1/2 , and ∆j+1/2 is defined as ∆j+1/2 = Qj+1 − Qj . δ is a small positive parameter. In this paper, the secondorder upwind TVD scheme suggested by Harten and Yee [10] is introduced, although there are a significant variety of TVD
3.1. Wave dispersion relations of plasma In order to confirm that the electromagnetic field is solved appropriately, we examine the wave dispersion relations of a cold plasma with periodic boundaries and a uniformly applied magnetic field. The results are compared to the analytic solutions. Linear perturbation theory reveals three waves of a cold plasma under the assumption of Darwin approximation. The R wave (electron cyclotron wave or Whistler wave) and the L wave (ion cyclotron wave) are circularly polarized transverse waves propagating parallel to the applied field, while the magnetoacoustic (MA) wave is a combination of longitudinal and transverse displacements with a propagation perpendicular to the applied field. Table 1 shows the non-dimensional parameters used for the simulation of the wave dispersion relation. The computational domain is one dimensional in the x direction. We initialize the uniform plasma and uniform applied field (Bx for the R and L waves and By for the MA wave) and observe the oscillation with time of the various field quantities. The distributions of the field quantities obtained here are stored and used for two-dimensional Fourier analyses in both time and space. The results are presented in Fig. 2: (a) analytic results, (b) numerical results for the R and L waves, and (c) numerical results for the MA wave. The analytic and numerical results are in good agreement. This dispersion relation does not depend obviously on the simulation code with or without the TVD scheme because the effect of the TVD scheme makes a difference only for a discontinuous surface of a magnetic field. 3.2. One-dimensional test simulation In order to estimate the characteristics of the computation method with the TVD scheme suggested in this paper, we carry out one-dimensional test simulations. Note that we divide the magnetic field into two parts: the induced plasma field (Bp ) and applied external field (Bd ). Fig. 3 shows the computational domain and working conditions. The one-dimensional, uniform, plasma flow is in the x direction from left to right. The plasma flow velocity is set to a value higher than the Whistler wave velocity. The applied field in the y direction, which is perpendicular to the plasma flow, is initialized at X = 40–50. The periodic boundary conditions are set at X = 0 and X = 100. This working condition poses serious obstacles to stably maintaining the calculation for hybrid PIC simulations. However, we choose such a serious condition
2030
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
(a) Theoretical.
(b) Numerical (parallel to B field).
(c) Numerical (perpendicular to B field).
Fig. 2. Wave dispersion relation for cold plasmas: (a) analytic results, (b) numerical results parallel to magnetic field, and (c) numerical results perpendicular to magnetic field.
Fig. 3. Computational domain and working conditions for one-dimensional test simulation.
on purpose to provide evidence of the effectiveness of the TVD scheme. Fig. 4 shows the magnetic field distributions at t = 2.5 when four kinds of hybrid PIC simulation code are used: (a) Terasawa’s scheme without smoothing, (b) the CAM-CL scheme without smoothing, (c) Terasawa’s scheme with smoothing as a primitive digital filter, and (d) Terasawa’s scheme with the TVD scheme. Here, the digital filter used in Fig. 4(c) is a simple binomial filter [13], as shown below. Bfiltered =
1 4
(Bj+1 + 2Bj + Bj−1 ).
(17)
The right-hand side of Eq. (17) is the output of Eq. (7), and the above filter function is adopted at every time step of the calculations. The solid line indicates the induced field, Bpy , and the dotted line indicates the applied magnetic field, Bdy , in Fig. 4. The value of β for the plasma, which represents the ratio of the plasma pressure to the magnetic pressure, is set to 1.0. Note that it is difficult to know the exact solution in the present one-dimensional test simulation because the governing equation system used in the hybrid PIC code is nonlinear. The plasma flow changes as a result of the interaction with the magnetic field. That is, if the plasma flow does not change by the interaction (or, without the interaction), the magnetic field is just advected with preservation of the waveform. In Fig. 4(a) and (b), an induced field that is associated with the diamagnetic current is generated and cancels out the applied field
at X = 40–50. Alternatively, the applied field is almost frozen into plasmas and convects downstream, as shown at about X = 60–80. Furthermore, short-wavelength oscillations appear at X = 20–60. These oscillations have a purely numerical origin. In fact, these are numerical non-physical oscillations originating from the evaluation method when the spatial derivative of the convection term of the induction equation is discretized numerically. The discretization method used in Fig. 4(a) and (b) is a second-order central difference scheme, which does not include any smoothing methods for the magnetic field. Terasawa’s scheme and the CAM-CL scheme produce almost the same results, because the difference between them involves the method used to discretize the time derivative term, rather than the spatial derivative term, of the magnetic induction equation. Numerical oscillations originating from the discretization of the spatial derivative are produced in both schemes. Such numerical oscillations may cause a decrease in the computational accuracy, as well as numerical instability and calculation failure. In Fig. 4(c), there are no oscillations in the numerical results, even though the waveform of the magnetic field is uniformly attenuated, because of the smoothing of the digital filtering. In general, such a filtering should be controlled by some empirical parameters and adopted locally and temporarily. The result shown in Fig. 4(c) is obviously overdamped. Thus, it may not be an appropriate solution. In Fig. 4(d), on the other hand, the numerical oscillations that would ordinarily be produced in the second-order central difference solution are completely prevented by the TVD scheme without overdamping the waveform of the magnetic field. In this case, the computation time increases by about 10–20% compared to the code without the TVD scheme. Note that the addition of the primitive smoothing shown in Fig. 4(c) also suppresses these oscillations but does not totally eliminate them. In this sense, filtering acts after the oscillations are produced by the basic numerical scheme. Fig. 5 shows the magnetic field distributions at t = 52.5 when the three kinds of hybrid PIC simulation code are used: the green line is Terasawa’s scheme without smoothing, the blue line is Terasawa’s scheme with the smoothing as a primitive digital filter, and the red line is Terasawa’s scheme with the TVD scheme. In the case without smoothing, numerical oscillations appear across the entire computational domain. Furthermore, in the case with smoothing as a digital filter, the waveform of the induced magnetic field becomes flatly broadened. Note that such a smoothing should be
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
2031
Fig. 4. Applied and induced magnetic field distributions (Bdy and Bpy ) at t = 2.5 for one-dimensional test simulation. (a) Terasawa’s scheme without smoothing, (b) CAM-CL scheme without smoothing, (c) Terasawa’s scheme with smoothing by a digital filter, and (d) Terasawa’s scheme with the TVD scheme.
Fig. 5. Magnetic field distributions (Bdy and Bpy ) at t = 52.5.
essentially applied locally and temporarily by some empirical parameters, as mentioned above. In this case, the magnetic field may produce overdamping because the digital filter used here is a simple binomial filter and is adopted at every time step. On the other hand, in the case with the TVD scheme, the magnetic field convects while keeping the waveform, which is slightly attenuated. Fig. 6 shows the temporal evolution of the magnetic energy of the system. Here, the magnetic energy of the system, EB , is defined as follows. EB =
(B2xj + B2yj + B2zj )1xj ,
(18)
j
where subscript j indicates the grid number. In Fig. 6, all of the
Fig. 6. Temporal evolution of magnetic energy of the system.
simulations show a decrease in the magnetic energy over time. In particular, the magnetic energy of the case with smoothing is greatly reduced compared to the other two cases. In fact, the decrease in magnetic energy corresponds to the attenuation of the waveform of the magnetic field. In the present simulations, the second term (electron pressure gradient term) and third term (Hall term) on the right-hand side of Eq. (7) are negligible compared to the first term (convection term). Here, if the plasma flow velocity is assumed to be uniform, the induction equation has almost the same form (the hyperbolic type) as shown below.
2032
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
Fig. 7. Space–time (x–t) diagram of induced magnetic field solved by Terasawa’s scheme with the TVD scheme.
∂ By ∂ By + Vx = 0. ∂t ∂x
(19)
Actually, the plasma flow changes as a result of the interaction with the magnetic field. That is, if the plasma flow does not change by the interaction (or, when the interaction is so weak), the magnetic field is just advected with preservation of the waveform and no magnetic energy decays. Note that the waveform of the magnetic field is absolutely oscillated or attenuated when a hyperbolic equation such as Eq. (19) is solved numerically.
In the case without smoothing, although the magnetic energy seems to be well conserved in Fig. 6, this is incorrect. In this case, the numerical oscillation of the magnetic field appears as shown in Fig. 4. Note that the conservation of the magnetic energy in the case without smoothing includes the influence of the numerical oscillation of the magnetic field. On the other hand, in the case with the TVD scheme, the waveform of the magnetic field is slightly attenuated over time because the TVD scheme also obviously includes the numerical dissipation. However, the decay time is longer than that with the digital filter. Fig. 7 shows the space–time (x–t) diagram of the induced field when Terasawa’s scheme with the TVD scheme is used. The convection of the frozen magnetic field can be observed over a long period, and the hybrid PIC code with the TVD scheme stabilizes the simulation numerically. This result provides evidence of the effectiveness and robustness of the TVD scheme. A parametric study should be done to confirm the validation of the hybrid PIC code with the TVD scheme. Fig. 8(a) and (b) show the time evolution of the induced field at t = 2.5, 5.0, and 7.5 under plasma conditions β = 0.01 and 10.0, respectively. Here, as above, β indicates the ratio of the plasma pressure to the magnetic pressure, and the applied field is set at X = 10–20. The induced field when β = 0.01 convects downstream, while it becomes diffuse compared with that when β = 10. In other words, the diffusion rate of the magnetic field increases with decreasing β . This is because of the difference in the thermal velocity of the ion particles. Fig. 9(a) and (b) show the ion flow velocity, ion number density, and induced field distribution at t = 5.0 when β = 0.01 and 10.0, respectively. These figures correspond to Fig. 8(a) and (b) at t = 5.0. In Fig. 9(a), the ion flow velocity and ion number
Fig. 8. Applied and induced magnetic field distributions for one-dimensional test simulations under conditions of (a) β = 0.01 and (b) β = 10.0.
Fig. 9. Ion number density, ion fluid velocity, and induced magnetic field distributions for one-dimensional test simulations under conditions (a) β = 0.01 and (b) β = 10.0.
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
2033
Table 2 Working conditions for two-dimensional simulation.
Fig. 10. Computational domain and working conditions for the two-dimensional simulation.
density change in the region X = 55–75. In the case of a low β (β = 0.01), the change in the plasma flow dynamics becomes sensitive to the interaction with the magnetic field because the thermal velocity is nearly zero. The induced field absolutely follows the plasma flow dynamics. Therefore, diffusion of the induced field occurs. In contrast, in the case of a high β (β = 10.0), the change in the plasma flow caused by the magnetic field vanishes into the noise of the thermal velocity, and the convection velocity of the induced field becomes constant and uniform. 3.3. Two-dimensional simulation In order to provide a simulation that is suitable for practical use, the global structure of a collision-less bow shock is examined by using the two-dimensional hybrid PIC code with the TVD scheme. Fig. 10 shows the computational domain. Solar-wind plasma is continuously injected from the left boundary, and the three remaining external boundaries are open to the flow of the plasma and electromagnetic fields. The applied current that creates the magnetic line dipole field is placed within the center of the computational domain. All of the particles reaching the outside of the simulation domain are lost from the system. All three components of the plasma velocity and electromagnetic fields are retained, but the computational domain is two dimensional in space, with d/dz = 0. Table 2 shows the working conditions. Here, rL /L is defined as the ratio of the ion Larmor radius of the solar wind at the magnetopause, rL , to a representative length of the magnetosphere, L. Fig. 11 shows the distributions of the ion number density and magnetic field lines at a steady state
Quantity and symbol
Value
rL /L Size of magnetosphere, L Coil current for line dipole Ion Larmor radius at the magnetopause, rL Plasma component Flow velocity, Vi Number density, Nc Ion and electron temperature, Ti , Te IMF strength Grid step size, dx Number of particles per unit cell Grid number (X –Y )
0.1 1000 km 390 kA 100 km proton + electron 400 km/s 5 × 106 l/m3 10 eV 3 nT 1 c /ωpi (=100 km) 36 300–300
by (a) Terasawa’s scheme without smoothing and (b) Terasawa’s scheme with the TVD scheme. In Fig. 11(a), numerical oscillation originating from the discretization of the magnetic induction equation appears in the entire region, and the structure of the bow shock is blurry. On the other hand, in Fig. 11(b), the structure of the bow shock and magnetic field lines is clearly captured, although no filtering is explicitly adopted in the simulation, except for the TVD scheme. 4. Summary In order to improve the numerical stability and robustness, a discretization procedure for the TVD scheme was introduced to the convection term of a magnetic induction equation in an electromagnetic hybrid particle-in-cell (PIC) plasma simulation code. The non-physical numerical oscillations that would ordinarily be produced in a solution estimated by a second-order central difference scheme at the discontinuities of magnetic fields were drastically prevented by the TVD scheme in the numerical results of one-dimensional test simulations, although the computation time increased by about 10–20% compared to the code without the TVD scheme. Furthermore, a two-dimensional simulation of the global structure of a collision-free bow shock, which is suitable for practical use, showed that the bow shock was clearly captured by using the hybrid PIC code with the TVD scheme. Acknowledgments This work was supported by the Japan Science and Technology Agency, Core Research for Evolutional Science and Technology (JST-CREST) under the project ‘‘High Performance Computing for Multi-Scale and Multi-Physics Phenomena’’.
Fig. 11. Distributions of ion number density and magnetic field lines of the two-dimensional simulation: (a) without smoothing, (b) with the TVD scheme.
2034
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
References [1] D. Winske, L. Yin, N. Omidi, H. Karimabadi, K. Quest, in: J. Büchner, C.T. Dum, M. Scholer (Eds.), Space Plasma Simulation, Springer-Verlag, Berlin, 2003, p. p. 136. [2] A.S. Lipatov, The Hybrid Multiscale Simulation Technology. An Introduction with Application to Astrophysical and Laboratory Plasmas, Springer-Verlag, Berlin, 2002. [3] D.S. Harned, J. Comput. Phys. 47 (1982) 452. [4] T. Terasawa, M. Hoshino, J. Sakai, T. Haba, J. Geophys. Res. 91 (A4) (1986) 4171. [5] D. Winske, K.B. Quest, J. Geophys. Res. 93 (A9) (1988) 9681. [6] E.J. Horowitz, D.E. Schumaker, D.V. Anderson, J. Comput. Phys. 84 (1989) 279.
[7] A.P. Matthews, J. Comput. Phys. 112 (1994) 102. [8] L. Gargaté, R. Bingham, R.A. Fonseca, L.O. Silva, Comput. Phys. Comm. 176 (2007) 419. [9] J. Müller, S. Simon, U. Motschmann, J. Schüle, K.H. Glassmeier, G.J. Pringle, Comput. Phys. Comm. 182 (2011) 946. [10] H.C. Yee, A Class of High-Resolution Explicit and Implicit Shock-Capturing Method, NASA Technical Memorandum 101080. [11] C.B. Laney, Computational Gasdynamics, Cambridge University Press, Cambridge, 1998. [12] A. Harten, J. Comput. Phys. 49 (1983) 357. [13] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, Taylor & Francis, 2005.
Contents lists available at SciVerse ScienceDirect
Computer Physics Communications journal homepage: www.elsevier.com/locate/cpc
Application of a total variation diminishing scheme to electromagnetic hybrid particle-in-cell plasma simulation Masaharu Matsumoto a,∗ , Yoshihiro Kajimura b , Hideyuki Usui a , Ikkoh Funaki b , Iku Shinohara b a
Department of Computational Science, Kobe University, Japan
b
Institute of Space and Astronautical Science, Japan Aerospace Exploration Agency, Sagamihara, Japan
article
info
Article history: Received 18 October 2011 Received in revised form 13 April 2012 Accepted 24 April 2012 Available online 3 May 2012 Keywords: Total variation diminishing Particle-in-cell code Electromagnetic hybrid plasma code Space plasma physics
abstract A discretization procedure for a total variation diminishing (TVD) scheme is introduced to an electromagnetic hybrid particle-in-cell (PIC) plasma simulation code in order to improve the numerical stability and resolution when calculating the plasma flow field in which magnetic field discontinuities (for example, Rankine–Hugoniot jump conditions for shock waves) are generated. In the hybrid PIC code used in this study, ions are treated as particles and electrons are assumed to be an inertia-less (massless) fluid. In the numerical results of one-dimensional test simulations, the TVD scheme significantly prevents non-physical, numerical oscillations, which would ordinarily be produced in the solution when the convection term of the magnetic induction equation in the hybrid PIC code is discretized by central difference schemes at magnetic field discontinuities. Furthermore, a two-dimensional simulation of the global structure of a collision-less bow shock, which is suitable for practical use, makes it possible to clearly capture the bow shock by using the hybrid PIC code with the TVD scheme. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Hybrid particle-in-cell (PIC) plasma simulation models [1,2], in which the various components of a plasma are treated in different ways and the plasma can be coupled to the electromagnetic field, have been widely utilized in plasma science and engineering over the past three decades. In the most common type of hybrid PIC code, the ions are treated kinetically as particles, the electrons are modeled as an inertia-less (mass-less) fluid, and the electromagnetic field is applied to the low-frequency (Darwin) approximation. Such models are designed to reveal the plasma phenomena that occur on shorter time and space scales than can be described by magnetohydrodynamics (MHD) models, although the detailed physical processes on electron inertial scales (for example, Debye length and electron Rarmar radius scales, and electron plasma frequency and electron cyclotron frequency time scales) are not resolved. In other words, hybrid PIC codes focus mainly on the plasma behaviors that occur on ion inertial scales (for example, ion Rarmar radius scales and ion cyclotron frequency time scales), namely, between the MHD and electron inertial scales. These models are used to reasonably simulate plasma phenomena
∗ Correspondence to: Department of Computational Science, Graduate School of System Infomatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan. Tel.: +81 0 78 803 6210; fax: +81 0 78 803 6210. E-mail addresses: [email protected] (M. Matsumoto), [email protected] (H. Usui). 0010-4655/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2012.04.021
in which ion kinetic effects are important. In hybrid PIC codes, the momentum equations for ion particles, generalized Ohm’s law (or momentum equation for an electron fluid), energy equation for an electron fluid, and Faraday’s induction equation (or induction equation for a magnetic field) may be solved numerically. Numerical algorithms and implementations of hybrid PIC codes have been developed and may be divided into several methods: the predictor–corrector scheme suggested by Harned [3], the model suggested by Terasawa et al. [4], the moment method suggested by Winske et al. [5], the QN3D method suggested by Horowitz et al. [6], the current advance method and cyclic leapfrog (CAM-CL) suggested by Matthews [7], amongst others. The predictor–corrector scheme [3] continues to be commonly utilized, even though it was first proposed as a hybrid PIC code. This is because it provides good energy conservation and robustness. However, it is well known that a considerable amount of shortwavelength Whistler noise is generated in this scheme. In order to avoid such noise, subcycling time steps are introduced and different time increments are used for the particles and fields in Terasawa’s model [4]. In other words, the time increment used for advancing the field, which is smaller than that for advancing the particles, is set to a value that is small enough to satisfy the requirement of the characteristic velocity of the Whistler wave, and the field calculations iterate until the particles are advanced. A second-order rational Runge–Kutta method is used to advance both the fields and the particles in this method. The moment method [5] has been proposed as a variation of Terasawa’s model.
2028
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
There are also subcycling time steps, but with a fourth-order Runge–Kutta method. Here, the ion fluid velocity is advanced using an MHD momentum equation. Instead of using the MHD momentum equation for advancing the ion fluid velocity, a fourthorder Adams–Bashforth extrapolation is also used [1]. In the QN3D method [6], functional iteration for advancing the fields is introduced instead of the subcycling time steps. In the CAM-CL scheme [7] that has becomes popular in recent years, the ion fluid velocity is calculated by advancing an extra half time step using a mixed level evaluation of the electric field. Additionally, with the rapid advancement of computer power, hybrid PIC codes have recently been developed that are suitable for massively parallel computation [8] or adaptive mesh refinement [9]. Although many variations of hybrid PIC codes obviously exist, with the exception of the methods mentioned above, the most notable differences between these codes involve the method used to advance the time steps or the methods used for the time integration of the particles and fields. In contrast, insufficient attention has been paid to the discretization of the spatial derivatives that appear in the magnetic induction equation, compared to that of the time derivative. The magnetic induction equation is a type of nonlinear hyperbolic partial differential equation, which has the convective derivative term shown below.
∂ b ∂ f (b) + = 0, ∂t ∂x
(1)
where b and f indicate an arbitrary physical quantity and function, respectively. Unfortunately, as is well known, the numerical analysis of such a hyperbolic conservation equation is not amenable to straightforward time and space integration and may give rise to numerical instabilities and non-physical oscillations if it is integrated in a naive way. In particular, almost all of the standard spatial discretization methods yield numerical results with excessive non-physical oscillations when a discontinuous surface appears in the physical quantity, for example, the Rankine–Hugoniot jump conditions of a shock wave. Therefore, this type of equation requires special care to be integrated numerically. One way to overcome this problem is to introduce an artificial viscosity as a smoothing filter. The idea that stable computations can be accomplished through an artificial viscosity or dissipater was used in an early scheme of computational fluid dynamics (CFD). The artificial viscosity smoothes out the strong gradients and attenuates the spurious oscillation in the solution. However, the difficulty with the artificial viscosity approach is that it is hard to determine an appropriate form that introduces just enough dissipation without causing unnecessary smearing. Over the past two decades, in CFD (including MHD), total variation diminishing (TVD) schemes [10,11], which were first introduced by Harten [12], have been widely used. The TVD schemes seldom generate spurious oscillation for nonlinear hyperbolic conservation laws because these guarantee that the total variation of the numerical solution will not increase under the numerical process in time. Any numerical scheme that gives rise to such oscillations does not satisfy the TVD conditions. For example, the central difference schemes that have been widely used in the convection terms of induction equations in conventional hybrid PIC codes are not TVD schemes. Based on the above, the purpose of this paper is to improve the numerical stability of a hybrid PIC code by applying a TVD scheme to the induction equation. 2. Numerical implementations The intended plasma in this paper is a singly ionized, fully ionized, and collision-less plasma. Fig. 1 shows a timechart of the hybrid PIC code used in this study, which is almost the same as Terasawa’s scheme [4]. Here, X# , V# , N , V, E, and B
Fig. 1. Timechart for the algorithm of the hybrid PIC code used in this paper.
indicate the position of the ion particles, velocity of the particles, number density, flow velocity, electric field, and magnetic field, respectively. Subscript i indicates the physical value of the ions. The general layout of the algorithm for the hybrid PIC code used in this study is as follows. Because the details of this hybrid PIC code can be found in Refs. [1,2], only a brief description is given below. We begin with the positions of the ion particles at time n, their velocities at time n + 1/2, and the electric and magnetic fields at time n. In step 1, the positions of the ion particles are updated from time n to n + 1/2 and then to n + 1. In step 2, the ion density and fluid velocity are calculated from the positions and velocities of the ion particles using the PIC method. In step 3, the magnetic field is updated by solving a magnetic induction equation while the ion density and ion fluid velocity are fixed in time, and then the electric field is calculated from the updated magnetic field. Here, subcycling time steps are used to update the magnetic field. Finally, in step 4, the velocities of the ion particles are updated using the Boris algorithm. Note that, in this paper, we focus exclusively on solving the induction equation in the hybrid PIC code. Therefore, the algorithm of step 3 will mainly be discussed later. For details about the numerical algorithms used in the hybrid PIC code, please refer to Refs. [1,2]. If the Darwin approximation and quasi-neutrality can be adopted for the intended plasma, the momentum equation for a mass-less electron fluid (generalized Ohm’s law), Maxwell’s equations, and the charge neutrality equation are shown as follows. E + Ve × B +
∇ Pe 2Ne
∂B = −∇ × E ∂t j=∇ ×B ∇ ·B=0 N = Ni = Ne ,
=0
(2) (3) (4) (5) (6)
where N , B, V, E, j, and Pe are the number density, magnetic field, flow velocity, electric field, current density, and electron pressure, respectively. Subscripts e represent electrons. Furthermore, the following normalizations are used: the number density, N, is normalized by the initial uniform number density, Nc ; the magnetic field, B, is normalized by the constant field magnitude, Bc ; the velocity, V, is normalized by the characteristic Alfven velocity, VA (=Bc /[µ0 mi Nc ]1/2 , where µ0 and mi are the vacuum permeability and ion mass, respectively); the electric field, E, is normalized by Bc VA ; the current, j, is normalized by eNc VA (e is the electric charge); and the pressure, P, is normalized by the magnetic pressure, B2c /2µ0 . The units for time and space are the reciprocal
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
of the ion cyclotron frequency, ωci−1 , and the ion inertial length, VA /ωci , respectively. The magnetic induction equation is given by substituting Eq. (2) into Eq. (3):
∂B ∇ Pe 1 = ∇ × Vi × B + − (∇ × B) × B . ∂t 2N N
∂Q ∂F + =0 ∂ t ∂ x By Vx By − Vy Bx Q= , F= . Bz Vx Bz − Vz Bx
(7)
1t n (F˜ − F˜jn−1/2 ). 1x j+1/2
(8)
(9)
Let Qjn be the numerical approximation of By or Bz at x = j1x and t = n1t, where 1x is the spatial mesh size and 1t is the time step. F˜ is a numerical flux vector function. The numerical flux for the second-order upwind TVD scheme suggested by Harten and Yee can be written as F˜j+1/2 =
1 2
(Fj+1 + Fj + φj+1/2 )
(10)
φj+1/2 = σ (cj+1/2 ) (gj + gj+1 ) − ψ(cj+1/2 + γj+1/2 )∆j+1/2
(11)
gj = minmod (∆j+1/2 , ∆j−1/2 )
(12)
minmod (x, y) = sgn(x) · max{0, min[|x|, y sgn(x)]}
(13)
γj+1/2
gj+1 − gj = σ (cj+1/2 ) ∆ 0 j+1/2
σ (z ) = ψ(z ) =
1 2
1t 2 ψ(z ) − z 1x
|z | (z 2 + δ 2 )/2δ
∆j+1/2 ̸= 0
(14)
∆j+1/2 = 0
|z | ≥ δ |z | < δ.
Quantity
Value
Time step Grid size Particles Grid number Ion density Magnetic field
0.01 0.5 36/cell 1024 1.0 1.0
schemes. Because the details of TVD schemes can be found in Ref. [10], only a brief description is given in this paper. Extensions to multi-dimensional equations can be obtained using a method similar to that mentioned above. The remaining second and third terms on the right-hand side of the magnetic induction equation (Eq. (7)) are discretized by a second-order central difference scheme. Subcycling time steps are also used [1,2,4]. That is, different time increments are used for update of the particles and fields. 3. Results
Here, note that the equation for Bx is ignored. We will discretize Eq. (8) while the ion fluid velocities (Vx , Vy , Vz ) are fixed in time. Generally, the numerical solvers for a nonlinear system equation of the hyperbolic conservation law are complicated, because there are characteristic velocities based on the number of equations. However, the system equation of Eq. (8) can be treated as two independent scalar equations for By and Bz , because the characteristic velocities of Eq. (8) are multiple roots of Vx . Many existing conservative schemes for the system of Eq. (8) use forward Euler time discretization. For simplicity, this time discretization is used for the discussion in this paper, and the scheme can be written as Qjn+1 = Qjn −
Table 1 Non-dimensional parameters for wave dispersion relation analysis.
Here, Eqs. (4) and (6) and the definition for the current, j = Ni Vi − Ne Ve , are substituted into Eq. (2). The first term on the right-hand side in Eq. (7) represents the convection term, while the second and third terms represent the electron pressure gradient term and Hall term, respectively. Note that the sum of the left-hand side and the first term on the right-hand side is the same hyperbolic form as Eq. (1). We attempt to apply the TVD scheme to the convection term of the magnetic induction equation. For simplicity, a onedimensional induction equation in which the time derivative and convection terms are extracted is considered. Then, the nonlinear system equation for the magnetic induction equations can be obtained, as shown below.
2029
(15)
(16)
Here, cj+1/2 indicates the characteristic velocity. In this case, cj+1/2 is equal to Vxj+1/2 , and ∆j+1/2 is defined as ∆j+1/2 = Qj+1 − Qj . δ is a small positive parameter. In this paper, the secondorder upwind TVD scheme suggested by Harten and Yee [10] is introduced, although there are a significant variety of TVD
3.1. Wave dispersion relations of plasma In order to confirm that the electromagnetic field is solved appropriately, we examine the wave dispersion relations of a cold plasma with periodic boundaries and a uniformly applied magnetic field. The results are compared to the analytic solutions. Linear perturbation theory reveals three waves of a cold plasma under the assumption of Darwin approximation. The R wave (electron cyclotron wave or Whistler wave) and the L wave (ion cyclotron wave) are circularly polarized transverse waves propagating parallel to the applied field, while the magnetoacoustic (MA) wave is a combination of longitudinal and transverse displacements with a propagation perpendicular to the applied field. Table 1 shows the non-dimensional parameters used for the simulation of the wave dispersion relation. The computational domain is one dimensional in the x direction. We initialize the uniform plasma and uniform applied field (Bx for the R and L waves and By for the MA wave) and observe the oscillation with time of the various field quantities. The distributions of the field quantities obtained here are stored and used for two-dimensional Fourier analyses in both time and space. The results are presented in Fig. 2: (a) analytic results, (b) numerical results for the R and L waves, and (c) numerical results for the MA wave. The analytic and numerical results are in good agreement. This dispersion relation does not depend obviously on the simulation code with or without the TVD scheme because the effect of the TVD scheme makes a difference only for a discontinuous surface of a magnetic field. 3.2. One-dimensional test simulation In order to estimate the characteristics of the computation method with the TVD scheme suggested in this paper, we carry out one-dimensional test simulations. Note that we divide the magnetic field into two parts: the induced plasma field (Bp ) and applied external field (Bd ). Fig. 3 shows the computational domain and working conditions. The one-dimensional, uniform, plasma flow is in the x direction from left to right. The plasma flow velocity is set to a value higher than the Whistler wave velocity. The applied field in the y direction, which is perpendicular to the plasma flow, is initialized at X = 40–50. The periodic boundary conditions are set at X = 0 and X = 100. This working condition poses serious obstacles to stably maintaining the calculation for hybrid PIC simulations. However, we choose such a serious condition
2030
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
(a) Theoretical.
(b) Numerical (parallel to B field).
(c) Numerical (perpendicular to B field).
Fig. 2. Wave dispersion relation for cold plasmas: (a) analytic results, (b) numerical results parallel to magnetic field, and (c) numerical results perpendicular to magnetic field.
Fig. 3. Computational domain and working conditions for one-dimensional test simulation.
on purpose to provide evidence of the effectiveness of the TVD scheme. Fig. 4 shows the magnetic field distributions at t = 2.5 when four kinds of hybrid PIC simulation code are used: (a) Terasawa’s scheme without smoothing, (b) the CAM-CL scheme without smoothing, (c) Terasawa’s scheme with smoothing as a primitive digital filter, and (d) Terasawa’s scheme with the TVD scheme. Here, the digital filter used in Fig. 4(c) is a simple binomial filter [13], as shown below. Bfiltered =
1 4
(Bj+1 + 2Bj + Bj−1 ).
(17)
The right-hand side of Eq. (17) is the output of Eq. (7), and the above filter function is adopted at every time step of the calculations. The solid line indicates the induced field, Bpy , and the dotted line indicates the applied magnetic field, Bdy , in Fig. 4. The value of β for the plasma, which represents the ratio of the plasma pressure to the magnetic pressure, is set to 1.0. Note that it is difficult to know the exact solution in the present one-dimensional test simulation because the governing equation system used in the hybrid PIC code is nonlinear. The plasma flow changes as a result of the interaction with the magnetic field. That is, if the plasma flow does not change by the interaction (or, without the interaction), the magnetic field is just advected with preservation of the waveform. In Fig. 4(a) and (b), an induced field that is associated with the diamagnetic current is generated and cancels out the applied field
at X = 40–50. Alternatively, the applied field is almost frozen into plasmas and convects downstream, as shown at about X = 60–80. Furthermore, short-wavelength oscillations appear at X = 20–60. These oscillations have a purely numerical origin. In fact, these are numerical non-physical oscillations originating from the evaluation method when the spatial derivative of the convection term of the induction equation is discretized numerically. The discretization method used in Fig. 4(a) and (b) is a second-order central difference scheme, which does not include any smoothing methods for the magnetic field. Terasawa’s scheme and the CAM-CL scheme produce almost the same results, because the difference between them involves the method used to discretize the time derivative term, rather than the spatial derivative term, of the magnetic induction equation. Numerical oscillations originating from the discretization of the spatial derivative are produced in both schemes. Such numerical oscillations may cause a decrease in the computational accuracy, as well as numerical instability and calculation failure. In Fig. 4(c), there are no oscillations in the numerical results, even though the waveform of the magnetic field is uniformly attenuated, because of the smoothing of the digital filtering. In general, such a filtering should be controlled by some empirical parameters and adopted locally and temporarily. The result shown in Fig. 4(c) is obviously overdamped. Thus, it may not be an appropriate solution. In Fig. 4(d), on the other hand, the numerical oscillations that would ordinarily be produced in the second-order central difference solution are completely prevented by the TVD scheme without overdamping the waveform of the magnetic field. In this case, the computation time increases by about 10–20% compared to the code without the TVD scheme. Note that the addition of the primitive smoothing shown in Fig. 4(c) also suppresses these oscillations but does not totally eliminate them. In this sense, filtering acts after the oscillations are produced by the basic numerical scheme. Fig. 5 shows the magnetic field distributions at t = 52.5 when the three kinds of hybrid PIC simulation code are used: the green line is Terasawa’s scheme without smoothing, the blue line is Terasawa’s scheme with the smoothing as a primitive digital filter, and the red line is Terasawa’s scheme with the TVD scheme. In the case without smoothing, numerical oscillations appear across the entire computational domain. Furthermore, in the case with smoothing as a digital filter, the waveform of the induced magnetic field becomes flatly broadened. Note that such a smoothing should be
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
2031
Fig. 4. Applied and induced magnetic field distributions (Bdy and Bpy ) at t = 2.5 for one-dimensional test simulation. (a) Terasawa’s scheme without smoothing, (b) CAM-CL scheme without smoothing, (c) Terasawa’s scheme with smoothing by a digital filter, and (d) Terasawa’s scheme with the TVD scheme.
Fig. 5. Magnetic field distributions (Bdy and Bpy ) at t = 52.5.
essentially applied locally and temporarily by some empirical parameters, as mentioned above. In this case, the magnetic field may produce overdamping because the digital filter used here is a simple binomial filter and is adopted at every time step. On the other hand, in the case with the TVD scheme, the magnetic field convects while keeping the waveform, which is slightly attenuated. Fig. 6 shows the temporal evolution of the magnetic energy of the system. Here, the magnetic energy of the system, EB , is defined as follows. EB =
(B2xj + B2yj + B2zj )1xj ,
(18)
j
where subscript j indicates the grid number. In Fig. 6, all of the
Fig. 6. Temporal evolution of magnetic energy of the system.
simulations show a decrease in the magnetic energy over time. In particular, the magnetic energy of the case with smoothing is greatly reduced compared to the other two cases. In fact, the decrease in magnetic energy corresponds to the attenuation of the waveform of the magnetic field. In the present simulations, the second term (electron pressure gradient term) and third term (Hall term) on the right-hand side of Eq. (7) are negligible compared to the first term (convection term). Here, if the plasma flow velocity is assumed to be uniform, the induction equation has almost the same form (the hyperbolic type) as shown below.
2032
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
Fig. 7. Space–time (x–t) diagram of induced magnetic field solved by Terasawa’s scheme with the TVD scheme.
∂ By ∂ By + Vx = 0. ∂t ∂x
(19)
Actually, the plasma flow changes as a result of the interaction with the magnetic field. That is, if the plasma flow does not change by the interaction (or, when the interaction is so weak), the magnetic field is just advected with preservation of the waveform and no magnetic energy decays. Note that the waveform of the magnetic field is absolutely oscillated or attenuated when a hyperbolic equation such as Eq. (19) is solved numerically.
In the case without smoothing, although the magnetic energy seems to be well conserved in Fig. 6, this is incorrect. In this case, the numerical oscillation of the magnetic field appears as shown in Fig. 4. Note that the conservation of the magnetic energy in the case without smoothing includes the influence of the numerical oscillation of the magnetic field. On the other hand, in the case with the TVD scheme, the waveform of the magnetic field is slightly attenuated over time because the TVD scheme also obviously includes the numerical dissipation. However, the decay time is longer than that with the digital filter. Fig. 7 shows the space–time (x–t) diagram of the induced field when Terasawa’s scheme with the TVD scheme is used. The convection of the frozen magnetic field can be observed over a long period, and the hybrid PIC code with the TVD scheme stabilizes the simulation numerically. This result provides evidence of the effectiveness and robustness of the TVD scheme. A parametric study should be done to confirm the validation of the hybrid PIC code with the TVD scheme. Fig. 8(a) and (b) show the time evolution of the induced field at t = 2.5, 5.0, and 7.5 under plasma conditions β = 0.01 and 10.0, respectively. Here, as above, β indicates the ratio of the plasma pressure to the magnetic pressure, and the applied field is set at X = 10–20. The induced field when β = 0.01 convects downstream, while it becomes diffuse compared with that when β = 10. In other words, the diffusion rate of the magnetic field increases with decreasing β . This is because of the difference in the thermal velocity of the ion particles. Fig. 9(a) and (b) show the ion flow velocity, ion number density, and induced field distribution at t = 5.0 when β = 0.01 and 10.0, respectively. These figures correspond to Fig. 8(a) and (b) at t = 5.0. In Fig. 9(a), the ion flow velocity and ion number
Fig. 8. Applied and induced magnetic field distributions for one-dimensional test simulations under conditions of (a) β = 0.01 and (b) β = 10.0.
Fig. 9. Ion number density, ion fluid velocity, and induced magnetic field distributions for one-dimensional test simulations under conditions (a) β = 0.01 and (b) β = 10.0.
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
2033
Table 2 Working conditions for two-dimensional simulation.
Fig. 10. Computational domain and working conditions for the two-dimensional simulation.
density change in the region X = 55–75. In the case of a low β (β = 0.01), the change in the plasma flow dynamics becomes sensitive to the interaction with the magnetic field because the thermal velocity is nearly zero. The induced field absolutely follows the plasma flow dynamics. Therefore, diffusion of the induced field occurs. In contrast, in the case of a high β (β = 10.0), the change in the plasma flow caused by the magnetic field vanishes into the noise of the thermal velocity, and the convection velocity of the induced field becomes constant and uniform. 3.3. Two-dimensional simulation In order to provide a simulation that is suitable for practical use, the global structure of a collision-less bow shock is examined by using the two-dimensional hybrid PIC code with the TVD scheme. Fig. 10 shows the computational domain. Solar-wind plasma is continuously injected from the left boundary, and the three remaining external boundaries are open to the flow of the plasma and electromagnetic fields. The applied current that creates the magnetic line dipole field is placed within the center of the computational domain. All of the particles reaching the outside of the simulation domain are lost from the system. All three components of the plasma velocity and electromagnetic fields are retained, but the computational domain is two dimensional in space, with d/dz = 0. Table 2 shows the working conditions. Here, rL /L is defined as the ratio of the ion Larmor radius of the solar wind at the magnetopause, rL , to a representative length of the magnetosphere, L. Fig. 11 shows the distributions of the ion number density and magnetic field lines at a steady state
Quantity and symbol
Value
rL /L Size of magnetosphere, L Coil current for line dipole Ion Larmor radius at the magnetopause, rL Plasma component Flow velocity, Vi Number density, Nc Ion and electron temperature, Ti , Te IMF strength Grid step size, dx Number of particles per unit cell Grid number (X –Y )
0.1 1000 km 390 kA 100 km proton + electron 400 km/s 5 × 106 l/m3 10 eV 3 nT 1 c /ωpi (=100 km) 36 300–300
by (a) Terasawa’s scheme without smoothing and (b) Terasawa’s scheme with the TVD scheme. In Fig. 11(a), numerical oscillation originating from the discretization of the magnetic induction equation appears in the entire region, and the structure of the bow shock is blurry. On the other hand, in Fig. 11(b), the structure of the bow shock and magnetic field lines is clearly captured, although no filtering is explicitly adopted in the simulation, except for the TVD scheme. 4. Summary In order to improve the numerical stability and robustness, a discretization procedure for the TVD scheme was introduced to the convection term of a magnetic induction equation in an electromagnetic hybrid particle-in-cell (PIC) plasma simulation code. The non-physical numerical oscillations that would ordinarily be produced in a solution estimated by a second-order central difference scheme at the discontinuities of magnetic fields were drastically prevented by the TVD scheme in the numerical results of one-dimensional test simulations, although the computation time increased by about 10–20% compared to the code without the TVD scheme. Furthermore, a two-dimensional simulation of the global structure of a collision-free bow shock, which is suitable for practical use, showed that the bow shock was clearly captured by using the hybrid PIC code with the TVD scheme. Acknowledgments This work was supported by the Japan Science and Technology Agency, Core Research for Evolutional Science and Technology (JST-CREST) under the project ‘‘High Performance Computing for Multi-Scale and Multi-Physics Phenomena’’.
Fig. 11. Distributions of ion number density and magnetic field lines of the two-dimensional simulation: (a) without smoothing, (b) with the TVD scheme.
2034
M. Matsumoto et al. / Computer Physics Communications 183 (2012) 2027–2034
References [1] D. Winske, L. Yin, N. Omidi, H. Karimabadi, K. Quest, in: J. Büchner, C.T. Dum, M. Scholer (Eds.), Space Plasma Simulation, Springer-Verlag, Berlin, 2003, p. p. 136. [2] A.S. Lipatov, The Hybrid Multiscale Simulation Technology. An Introduction with Application to Astrophysical and Laboratory Plasmas, Springer-Verlag, Berlin, 2002. [3] D.S. Harned, J. Comput. Phys. 47 (1982) 452. [4] T. Terasawa, M. Hoshino, J. Sakai, T. Haba, J. Geophys. Res. 91 (A4) (1986) 4171. [5] D. Winske, K.B. Quest, J. Geophys. Res. 93 (A9) (1988) 9681. [6] E.J. Horowitz, D.E. Schumaker, D.V. Anderson, J. Comput. Phys. 84 (1989) 279.
[7] A.P. Matthews, J. Comput. Phys. 112 (1994) 102. [8] L. Gargaté, R. Bingham, R.A. Fonseca, L.O. Silva, Comput. Phys. Comm. 176 (2007) 419. [9] J. Müller, S. Simon, U. Motschmann, J. Schüle, K.H. Glassmeier, G.J. Pringle, Comput. Phys. Comm. 182 (2011) 946. [10] H.C. Yee, A Class of High-Resolution Explicit and Implicit Shock-Capturing Method, NASA Technical Memorandum 101080. [11] C.B. Laney, Computational Gasdynamics, Cambridge University Press, Cambridge, 1998. [12] A. Harten, J. Comput. Phys. 49 (1983) 357. [13] C.K. Birdsall, A.B. Langdon, Plasma Physics via Computer Simulation, Taylor & Francis, 2005.