Structure-preserving geometric particle-in-cell algorithm suppresses finite-grid instability -- Comment on "Finite grid instability and spectral fidelity of the electrostatic Particle-In-Cell algorithm'' by Huang et al
aa r X i v : . [ phy s i c s . p l a s m - ph ] A p r Structure-preserving geometric particle-in-cell algorithmsuppresses finite-grid instability – Comment on “Finite gridinstability and spectral fidelity of the electrostatic Particle-In-Cellalgorithm” by Huang et al.
Jianyuan Xiao and Hong Qin
2, 1, ∗ School of Physical Sciences, University of Scienceand Technology of China, Hefei, 230026, China Plasma Physics Laboratory, Princeton University, Princeton, NJ 08543, U.S.A
Abstract
A recent paper by Huang et al. [Computer Physics Communications 207, 123 (2016)] thoroughlyanalyzed the Finite Grid Instability (FGI) and spectral fidelity of standard Particle-In-Cell (PIC)methods. Numerical experiments were carried out to demonstrate the FGIs for two PIC meth-ods, the energy-conserving algorithm and the momentum-conserving algorithm. The paper alsosuggested that similar numerical experiments should be performed to test the newly developedStructure-Preserving Geometric (SPG)-PIC algorithm. In this comment, we supply the results ofthe suggested numerical experiments, which show that the SPG-PIC algorithm is able to suppressthe FGI. ∗ [email protected] L × × L = 33, ∆ x = 1,and ω p = 2 π/L . The time step is set to ∆ t = 0 . ω − p . The numbers of sampling pointsper grid for both electron and ions are 300, and the mass ratio and charge ratio betweenelectrons and ions are 1 : 3672 and − δx ( x p ) = LA cos (2 πM x p /L ) / (2 πM ), and their velocities are v x ( x ) = 0 . LAω p sin (2 πM x/L ) / (2 πM ), where A = 0 .
01 and M = 9. Initial electric field is E x ( x ) = − Lq e A cos (2 πM x/L ) / (2 πM ).The simulation first is performed to t = 220 ω − p . The resulting mode spectrum, finalvelocity distribution, energy and momentum evolution are plotted in Figs. 1 and 2, whichcorrespond to Fig. 3 and Fig. 4 of Ref. [1], respectively. These results show that even thoughmode alias still exists in the SPG-PIC algorithm, the FGI is suppressed. While mode aliaseffect, as an error of spatial discretization, is inevitable in any spatial grid, its existencedoes not necessarily imply unstable numerical eigenmodes will be exited. Note that whena unstable numerical eigenmode is excited, all components of the dynamics, especially thedominant ones, of the discrete system will grow exponentially. Unfortunately, such FGIsdo exist in the standard PIC methods. As demonstrated in Ref. [1], in about 30 plasmaoscillation periods (200 ω − p ), the total energy error for the MC-PIC algorithm exceeds 200% , and the total momentum error for the EC-PIC algorithm exceeds 70%. On the other hand,Fig. 2 shows that for the SPG-PIC algorithm, the total energy error is less than 1% , and the2otal momentum error is less than 0 . . The observed suppressing of FGI for the SPG-PIC algorithm can be attributed to thestructure-preserving nature of its spatial discretization. The charge deposition and fieldinterpolation are derived from a variational principle using the techniques of Whitney inter-polation forms [5, 10, 14] or finite element discrete exterior calculus [9, 13], which preservesthe discrete gauge symmetry and the discrete exterior calculus structure of the electromag-netic field. As a result, physical laws, such as the charge conservation and ∇ · B = 0,are satisfied exactly by the discrete system. This result is consistent with Ref. [1]’s conclu-sion that charge deposition and field interpolation can be optimally designed to suppress orreduce FGIs.Another feature of the SPG-PIC algorithm is the preserving of non-canonical symplec-tic structure for time-integration, which in general bounds simulation errors on conservedquantities for a very long time. We run the simulation for 500 longer to t = 10000 ω − p , andthe result is plotted in Fig. 3. Over this long simulation time, the total energy error andtotal momentum error are bounded by 1% and 2%, respectively.We finish this comment with two footnotes. First, the SPG-PIC algorithm used is forVlasov-Maxwell system in 3D configuration space. For the simplified geometry and simula-tion parameters of the present numerical experiments, the dominated modes of the discretesystem are longitudinal electrostatic modes. Secondly, for the PAS methods, symplectictime-integration can also be adopted. For example, Cary and Doxas [15, 16] first applieda canonical symplectic algorithm to simulate the particle-and-mode Hamiltonian models[17–22] for the Vlasov-Poisson system. ACKNOWLEDGMENTS
We thank Prof. John Cary and Prof. Dominique Escande for fruitful discussions. This re-search is supported by the National Key Research and Development Program (2016YFA0400600,2016YFA0400601, 2016YFA0400602 and 2018YFE0304100), the National Natural ScienceFoundation of China (NSFC-11775219, NSFC-11575186 and NSFC-11805273), China Post-doctoral Science Foundation (2017LH002), the Fundamental Research Funds for the Central3 mode number050100150200 ω p t − − − − − − − − − tω pe − − − − − − − − l og ( | ff t( E ) | mode 9mode 15mode 6 Figure 1. Mode spectrum simulated by the SPG-PIC algorithm. Full mode spectrum log( | fft( E ) | )as a function of time is shown in the left figure, and evolution of the mode amplitude for modenumber k = 9 , , x/ ∆ x . . . . . . v e , x / c tω pe . . . . . . . . E/E (0)
P/P (0)
Figure 2. Electron distribution at t = 200 ω − p (left) and the evolution of energy and momentum(right). Universities (WK2030040096) and the U.S. Department of Energy (DE-AC02-09CH11466). [1] C.-K. Huang, Y. Zeng, Y. Wang, M. D. Meyers, S. Yi, and B. J. Albright, Computer PhysicsCommunications , 123 (2016).[2] C. K. Birdsall and A. B. Langdon,
Plasma Physics via Computer Simulation (IOP Publishing,1991) p. 293. tω pe − − − − − − − − l og ( | ff t( E ) | mode 9mode 15mode 6 tω pe . . . . . . . . E/E (0)
P/P (0)
Figure 3. Long term evolution of spectrum (left), total energy and momentum (right) simulatedby the SPG-PIC algorithm.[3] H. R. Lewis, Journal of Computational Physics , 136 (1970).[4] E. Evstatiev and B. Shadwick, Journal of Computational Physics , 376 (2013).[5] J. Squire, H. Qin, and W. M. Tang, Physics of Plasmas , 084501 (2012).[6] J. Xiao, J. Liu, H. Qin, and Z. Yu, Phys. Plasmas , 102517 (2013).[7] J. Xiao, H. Qin, J. Liu, Y. He, R. Zhang, and Y. Sun, Physics of Plasmas , 112504 (2015).[8] Y. He, H. Qin, Y. Sun, J. Xiao, R. Zhang, and J. Liu, Physics of Plasmas , 124503 (2015).[9] Y. He, Y. Sun, H. Qin, and J. Liu, Physics of Plasmas , 092108 (2016).[10] J. Xiao, H. Qin, P. J. Morrison, J. Liu, Z. Yu, R. Zhang, and Y. He, Physics of Plasmas ,112107 (2016).[11] H. Qin, J. Liu, J. Xiao, R. Zhang, Y. He, Y. Wang, Y. Sun, J. W. Burby, L. Ellison, andY. Zhou, Nuclear Fusion , 014001 (2016).[12] J. Xiao, H. Qin, J. Liu, and R. Zhang, Physics of Plasmas , 062112 (2017).[13] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker, Journal of Plasma Physics (2017).[14] J. Xiao, H. Qin, and J. Liu, Plasma Science and Technology , 110501 (2018).[15] J. Cary and I. Doxas, Journal of Computational Physics , 98 (1993).[16] I. Doxas and J. R. Cary, Physics of Plasmas , 2508 (1997).[17] H. E. Mynick and A. N. Kaufman, Physics of Fluids , 653 (1978).[18] D. F. Escande, Physica Scripta T2A , 126 (1982).
19] D. F. Escande, in
Hamiltonian Dynamical Systems , edited by R. S. MacKay and J. D. Meiss(Adam Hilger, 1987) pp. 446–461.[20] D. F. Escande, in
Nonlinear World , Vol. 2, edited by V. Bar’yakhtar, V. M. Chernousenko,N. S. Erokhin, A. G. Sitenko, and V. E. Zakharov (World Scientific, Singapore, 1989) pp.817–836.[21] D. F. Escande, in
Large Scale Structures in Nonlinear Physics , edited by J. Fournier andP. Sulem (Springer, Berlin, 1991) pp. 73–104.[22] D. F. Escande, D. Bénisti, Y. Elskens, D. Zarzoso, and F. Doveil, Reviews of Modern PlasmaPhysics (2018).(2018).