Verification of long wavelength electromagnetic modes with a gyrokinetic-fluid hybrid model in the XGC code

TL;DR

This paper extends the XGC1 gyrokinetic PIC code to include a hybrid model with fluid electrons, enabling verification of long wavelength electromagnetic modes like shear Alfvén and tearing modes in different magnetic geometries.

Contribution

It introduces a gyrokinetic-fluid hybrid model in XGC1 for electromagnetic mode verification, combining gyrokinetic ions with fluid electrons.

Findings

Successful verification of shear Alfvén waves.

Verification of resistive tearing modes.

Application in cylindrical and toroidal geometries.

Abstract

As an alternative option to kinetic electrons, the gyrokinetic total-f particle-in-cell (PIC) code XGC1 has been extended to the MHD/fluid type electromagnetic regime by combining gyrokinetic PIC ions with massless drift-fluid electrons analogous to Chen and Parker, Physics of Plasmas 8, 441 (2001). Two representative long wavelength modes, shear Alfv\'en waves and resistive tearing modes, are verified in cylindrical and toroidal magnetic field geometries.

Equations

$$\frac{\partial \delta n_e}{\partial t}$$

$$+\delta {\mathbf{B}}_{\perp }\cdot \nabla \frac{j_0}{eB}-{\mathbf{v}}_E\cdot \nabla \left(n_0+\delta n_e\right)$$

$$-\frac{2n_0}{B^3}\left(\mathbf{B}\times \nabla B\right)\cdot \nabla \varphi$$

$$+\frac{2}{e{B}^3}\left(\mathbf{B}\times \nabla B\right)\cdot \nabla \delta P_e,$$

$$\frac{\partial A_{\parallel }}{\partial t}$$

$$E_{\parallel }$$

$$-\nabla \cdot \frac{{ϵ}_0\chi }{e}\nabla \varphi =\delta n_i-\delta n_e,$$

Full text

This article describes the verification of two important MHD/fluid type, long-wavelength, electromagnetic modes after the addition of an optional kinetic-fluid hybrid model to the gyrokinetic particle-in-cell code XGC1 Ku, Chang, and Diamond (2009). This work complements – as a more economical alternative – the fully implicit, fully kinetic electromagnetic formulation, that is also being developed for XGC1 Ku et al. (2016a).

The importance of MHD/fluid type electromagnetic modes in magnetically confined fusion devices, which operate regularly at moderate to high $\beta=2\mu_{0}P/B^{2}$ (the ratio of thermodynamic to magnetic pressure), is widely recognized. Examples are neoclassical tearing modes Connor et al. (1988), sawtooth oscillations von Goeler, Stodiek, and Sauthoff (1974), and edge localized modes Huysmans (2005). Gyrokinetic electromagnetic codes such as GYRO Candy and Waltz (2003); Candy (2005), GS2 Kotschenreuther, Rewoldt, and Tang (1995), GENE Pueschel, Kammerer, and Jenko (2008), GEM Parker, Chen, and Kim (2000); Chen and Parker (2001, 2003) have been available with increasing physics capability for more than a decade and have also been used to study those modes. However, their application to long wavelength MHD/fluid type instabilities has been difficult, especially for PIC codes, due to the so called “cancellation problem” Cummings (1995); Mishchenko, Hatzky, and Könies (2004). Recently, several methods to overcome the cancellation problem with kinetic electrons have been developed for particle codes: the control variate method Hatzky, Könies, and Mishchenko (2007), a special splitting of the vector potential Mishchenko et al. (2014a, b) (used e.g. by the EUTERPE code), and split-weight methods Lee et al. (2001); Chen and Parker (2003) (used in GEM, and being further developed in GTS Startsev et al. (2016)). The XGC1 code Ku, Chang, and Diamond (2009) recently demonstrated fully kinetic electromagnetic capability without cancellation problem Ku et al. (2016a) using a fully implicit electromagnetic scheme based on the work by Chen and Chacón Chen and Chacón (2014, 2015); Chacón and Chen (2016). These methods are computationally expensive for long wave length MHD/fluid type modes even without the cancellation problem. The cheapest way to study these modes is to use fluid electrons instead of electron particles.

Long wavelength electromagnetic physics in the global edge region have so far been studied with fluid and MHD codes (some of them with ad-hoc kinetic ion effect) such as BOUT++ Dudson et al. (2009); Liu et al. (2014), M3D Fu et al. (2006); Breslau, Park, and Jardin (2006), M3D-C1 Ferraro and Jardin (2009); Ren et al. (2015), and JOREK Czarny and Huysmans (2008); Huysmans et al. (2010), which neglect important effects that drive the plasma to a non-thermal equilibrium. Since kinetic ion effects on fluid/MHD modes as well as microturbulence are expected to be important in the plasma edge region, e.g. for the physics of edge localized modes (ELMs), kinetic ballooning modes (KBM) and others, we will improve the fluid and MHD approach by coupling gyrokinetic ions to the massless electron fluid hybrid model utilized in the GEM code Parker, Chen, and Kim (2000); Chen and Parker (2001). Although the fluid treatment of the electrons drops some important effects such as the trapped electron mode (TEM), it is still attractive because its implementation is rather straightforward without the cancellation issue, low $k_{\perp}$ fluid/MHD modes are important for ELM activity, and it is economical with computing time.

The fluid-kinetic hybrid version of the XGC1 code used for this report combines gyrokinetic ions in the $\delta$f formalism Ku, Chang, and Diamond (2009) (which, if done correctly, can be made identical to the total-f formalism Ku et al. (2016b)) with massless drift-fluid electrons Chen and Parker (2001); Chen et al. (2015). The electron density continuity equation is given by:

[TABLE]

where $\bm{B}$ is the axisymmetric background magnetic field, $\delta\bm{B}_{\perp}=\nabla A_{\parallel}\times\hat{\bm{b}}$ is the perturbed magnetic field, $\hat{\bm{b}}=\bm{B}/B$, and $A_{\parallel}$ is the component of the perturbed vector potential along the background magnetic field; $\mu_{0}$ is the vacuum permeability; $u_{\parallel,i}=\int\mathrm{d}^{3}vv_{\parallel}\delta f_{i}/n_{0}$ is the parallel ion fluid flow; $\delta f_{i}$ is the perturbed ion guiding center distribution function; $n_{e}=n_{0}+\delta n_{e}$ is the electron density; $j_{0}=\hat{\bm{b}}\cdot\nabla\times(\bm{B}/\mu_{0})$ is the equilibrium current density; $\delta P_{e}=e\delta n_{e}T_{0,e}$ is the perturbed iso-thermal electron pressure, and $\bm{v}_{E}=\frac{1}{B}\hat{\bm{b}}\times\nabla\phi$ is the $\bm{E}\times\bm{B}$ drift. We also used the relation $u_{\parallel,e}=(\nabla_{\perp}^{2}A_{\parallel})/(e\mu_{0}n_{0})+u_{\parallel,i}$. The time evolution of the perturbed vector potential is given by the definition of the electric field and Ohm’s law,

[TABLE]

Here, $\tilde{\bm{b}}=\hat{\bm{b}}+\delta\bm{B}_{\perp}/B$, $P_{0,e}=en_{0}T_{0,e}$ is the background electron pressure, and $j_{\parallel}=en_{0}(u_{\parallel,i}-u_{\parallel,e})$. Finally, the gyrokinetic Poisson equation in the long wave length limit is

[TABLE]

where $\chi=(\rho_{i}/\lambda_{D,i})^{2}=(c/v_{A})^{2}$ is the ion electric susceptibility, $\rho_{i}$ is the ion gyro radius, $\lambda_{D,i}$ is the ion Debye length, and $v_{A}=B/(\mu_{0}m_{i}n_{0})^{1/2}$ is the Alfvén speed. The ion density is $\delta n_{i}=\int\mathrm{d}^{3}v\langle\delta f_{i}\rangle_{\rho}$, where $\langle\dots\rangle_{\rho}$ indicates gyro-averaging.

The massless electron approximation is valid in the limit $v_{A}/v_{t,e}\rightarrow 0$ or $\beta_{e}m_{i}/m_{e}\gg 1$, where $\beta_{e}=2\mu_{0}P_{e}/B^{2}$.

Both explicit and implicit time integrators have been implemented. A second order Runge-Kutta (RK2) method has been utilized for the time integration of the combined particle-fluid system for many of the results discussed in this work. In the first step, $\delta n_{e}(t+\Delta t/2)$ and $A_{\parallel}(t+\Delta t/2)$ are evaluated using, $\phi(t)$, $\delta n_{i}(t)$ and $u_{\parallel,i}(t)$. Then the particles are pushed for a half time step to evaluate $\delta n_{i}(t+\Delta t/2)$ and $u_{\parallel,i}(t+\Delta t/2)$. In the second step, we evaluate $\phi(t+\Delta t/2)$ and then push the particles for a full time step to obtain $\phi(t+\Delta t)$, $\delta n_{i}(t+\Delta t)$ and $u_{\parallel,i}(t+\Delta t)$.

Implicit time stepping methods have been implemented using the PETSc TS framework Balay et al. (2016a, b, 1997) to overcome the restrictions on the time step of explicit methods. The particle terms $\delta n_{i}$ and $u_{\parallel,i}$ are treated as non-linear contributions to the system of electron fluid equations and are fully integrated into PETSc’s nonlinear solver residual, but only the electron fluid terms are included in the Jacobian. The Newton method is used to solve the non-linear equations, which requires one particle push per evaluation of the residual.

For the verification of shear Alfvén wave physics, we use a minimal system that supports this mode: linearized equations versions of (1)-(3) with the closure $E_{\parallel}=\eta_{e}j_{\parallel}$. In addition, we neglect the terms related to the curvature and $\nabla B$ drift in Eq. (1). It is straightforward to prove in cylindrical geometry that the dispersion relation of the resulting reduced system yields $\omega=[v_{A}2k_{\parallel}^{2}-4(\eta_{e}/(2\mu_{0}))^{4}k_{\perp}^{4}]^{1/2}+i\eta_{e}/(2\mu_{0})k_{\perp}^{2}$. The first verification test of the shear Alfvén dispersion relation was conducted in cylindrical geometry with concentric, circular flux-surfaces with minor radius $a=1$, constant safety factor $q=3$, $\beta_{e}=1.5\cdot 10^{-2}$, and $\eta_{e}=10^{-6}\,\Omega\mathrm{m}$. The simulation was initialized with a global perturbation of $A_{\parallel}$ centered around $r/a=0.67$ containing toroidal mode numbers $n=1\dots 4$ and poloidal mode numbers $m=0\dots 4$. With this large scale variation in the radial and poloidal direction, the low resistivity does not influence the real frequency much but still serves as a check for the resistive dissipation of the reduced shear Alfvén wave system (with $k_{\perp}\sim 1/a$). The time step for this simulation was $\Delta t=1.36\cdot 10^{-8}\,\mathrm{s}\approx 10^{-2}\tau_{A}$, where $\tau_{A}=R_{0}/v_{A}$. The total duration of the simulation is $1.36\cdot 10^{-3}\,\mathrm{s}\approx 1000\tau_{A}$.

Figure 1 shows the shear Alfvén spectrum obtained from this simulation. The parallel wave number was determined as $k_{\parallel}=\hat{\bm{b}}\cdot\bm{k}=(B_{P}/B)k_{\theta}+(B_{T}/B)k_{\varphi}$. The mode frequency is the median of the intensity for each value of $k_{\parallel}$ and the error bars indicate the decay length of the mode intensity around its median. The increasing width of the error bars at $k_{\parallel}>0.5$ indicates decreasing overall intensity due to the low toroidal and poloidal mode numbers used to initialize the simulation. The steps in the frequency spectrum are an artifact of the interpolation of the intensity from $(k_{\theta},k_{\varphi})$ space to a common $k_{\parallel}$ scale.

Similar tests in toroidal geometry have been performed in a slightly modified version of the standard cyclone geometry, with $R_{0}=1.7$, $a/R_{0}=0.358$, $B_{0}=1.9$, constant $q=2$, and $T_{0}=2$ keV. The density is varied between $1.875\cdot 10^{19}\,\mathrm{m}^{-3}$ and $6\cdot 10^{20}\,\mathrm{m}^{-3}$ to achieve values of $\beta_{e}$ between $0.4$% and $13.4$%. The time step is $\Delta t\approx 5\cdot 10^{-2}\tau_{A}$ and the total simulation time is $\approx 40\tau_{A}$. Figure 2 a) shows a poloidal wave number scan ($m=6\dots 10$) of the frequency of the $n=4$ shear Alfvén wave. The numerical frequencies agree very well with the (approximate) analytical result $\omega\propto(2\pi/L_{\parallel})(n-m/q)v_{A}$, where $L_{\parallel}$ is the parallel connection length for one poloidal circuit. The deviations are caused by the variation of the field line pitch along magnetic field lines. We find that $\omega\tau_{A}$ is independent of $\beta_{e}$ as expected because only the density $n_{0}$ was varied in this test.

Since the kinetic-fluid hybrid approach is especially useful for the simulation of low-n tearing modes, we benchmarked the $(m,n)=(2,1)$ tearing mode in cylindrical and toroidal geometry against the GEM code Chen et al. (2015) and M3D-K Cai and Fu (2012), respectively. We do not consider the effect of kinetic ions in this benchmark. The only term added to the electron fluid equations compared to the terms kept in the shear Alfvén case is the kink drive $\delta\bm{B}_{\perp}\cdot\nabla[j_{0}/(eB)]$ in Eq. (1) to be consistent with GEM’s eigenvalue solver Chen et al. (2015).

For the benchmark against the GEM code, we use the case described in Ref. Chen et al., 2015: concentric, circular flux-surfaces in cylindrical geometry, $R_{0}=1.7$ m, $a=0.425$ m ($R_{0}/a=4$), $B_{0}=1.906$ T, $q=1.5[1+(r/a)^{2}]$, $Z=1$, $m_{i}/m_{p}=2.5$, and constant density $n_{0}=3.886\cdot 10^{20}\,\mathrm{m}^{-3}$. Since the electron fluid equations used for this benchmark have no temperature dependence, we can use a constant temperature profile $T_{0,e}=45.63$ eV, which yields the same on-axis $\beta_{e}$ of $4\cdot 10^{-3}$ and relative domain size $a/\rho_{i}\approx 740$ as in Ref. Chen et al., 2015. The resonant surface for the $(2,1)$ tearing mode is at $(r/a)_{c}\approx 0.577$ corresponding to the normalized poloidal magnetic flux $\psi_{N,c}=0.411$. In order to be able to resolve the resonance layer of the $(2,1)$ tearing mode also at low resistivity, the radial resolution of our computational mesh varies between $0.5$ mm around the resonant surface to a maximum of $8$ mm far from the resonant surface. The relations between the normalized resistivity $\eta_{N}$ and growth rate $\gamma_{N}$ used in Ref. Chen et al., 2015 and the corresponding values $\eta_{e}$ and $\gamma$ in SI-units are $\eta_{N}=(en_{0}/B_{0})\eta_{e}$ and $\gamma_{N}=m_{p}/(eB_{0})\gamma$, where $m_{p}$ is the proton mass. The results of a resistivity scan of the growth rate of the $(2,1)$ tearing mode in this geometry is shown in Fig. 3. The growth rates evaluated with XGC1 show excellent agreement with the growth rates computed with GEM’s eigenvalue solver that uses the MHD approximation for the ion polarization density (Fig. 3 in Ref. Chen et al., 2015). We did not include an XGC1 data point for $\eta_{N}=10^{-7}$ because of the very strict resolution requirements of about $2.5\cdot 10^{-4}$ m or less for this low resistivity. Using the Crank-Nicolson method, the implicit time integrator could speed up these simulation by a factor of more than 10. For $\eta_{N}=10^{-6}$ a time step of $\Delta t=2.7\tau_{A}$ could be used.

For the benchmark against M3D-K in toroidal geometry, we use a Grad-Shafranov equilibrium generated with the FLOW code Guazzotto et al. (2004) with a fixed circular boundary, $R_{0}=5.76$ m, $a=1$ m, $B_{0}=1$ T, $q=1.5+2\psi_{N}^{2}$, $m_{i}/m_{p}=2.5$, and constant $n_{0}=10^{20}\,\mathrm{m}^{-3}$ and $T_{0,e}=100$ eV, so that $\beta_{e}=4\cdot 10^{-3}$ and $\beta_{e}m_{i}/m_{e}=18.4$. The resonant surface of the $(2,1)$ tearing mode is located at $\psi_{N}=0.5$. The radial resolution of the computational mesh is $1.5$ mm between approximately $\psi_{N}=0.4$ and $\psi_{N}=0.6$ and up to $1.2$ cm away from the tearing layer. For the normalized resistivity of $\eta_{M3D}=10^{-4}$ used in Ref. Cai and Fu, 2012 and the corresponding normalization relations, we obtain a resistivity in SI units of $\eta_{e}=\mu_{0}(a^{2}/\tau_{A})\eta_{M3D}=3.01\cdot 10^{-4}$. The XGC1 growth rate calculation used a time step of $\Delta t=7\cdot 10^{-3}\tau_{A}$ and ran for a total time of approximately $350\,\tau_{A}$. Figures 4 a)-d) show the mode structure of the growing $(2,1)$ mode, which exhibit the usual tearing structure. For comparison with Ref. Cai and Fu, 2012, we use reduced MHD quantities, the perturbed current $Ru_{\parallel,e}$, and the velocity stream function $\phi/B$.

The growth rate we obtain from the XGC1 calculation is $\gamma=1.12\cdot 10^{-2}\tau_{A}^{-1}$ and compares well to Ref. Cai and Fu, 2012. The relative difference between the XGC1 and the M3D-K result $6$%.

In order to add gyrokinetic ion effects to electromagnetic fluid/MHD instabilities, the gyrokinetic edge turbulence code XGC1 has been modified by replacing the kinetic electrons by massless drift-fluid electrons Parker, Chen, and Kim (2000); Chen and Parker (2001). Explicit and implicit time integration methods have been implemented and tested. We verified shear Alfvén wave physics against the analytical solution and benchmarked the massless fluid model for resistive tearing modes against the codes GEM and M3D-K. The hybrid model in XGC will be further developed into a total-f code with the aim of studying the onset of edge localized modes across the magnetic separatrix surface. Verification of the kinetic version of peeling-ballooning modes, and kinetic ballooning modes will be reported in a subsequent paper.

The authors would like to thank Guoyong Fu, Stephen Abbott, and Eduardo D’Azevedo for fruitful discussions. Support for this work was provided through the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy Office of Advanced Scientific Computing Research and the Office of Fusion Energy Sciences. Notice: This manuscript has been authored by Princeton University under Contract No. DE-AC02-09CH11466 with the U.S. Department of Energy. The publisher, by accepting the article for publication, acknowledges that the United States Government retains a non- exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

• Ku, Chang, and Diamond (2009) S. Ku, C. Chang,  and P. Diamond, Nucl. Fusion 49, 115021 (2009).
• Ku et al. (2016a) S.-H. Ku, R. Hager, C. Chang, L. Chacón,  and G. Chen, in B. Am. Phys. Soc. (American Physical Society, 2016) p. GP10.00112.
• Connor et al. (1988) J. W. Connor, S. C. Cowley, R. J. Hastie, T. C. Hender, A. Hood,  and T. J. Martin, The Physics of Fluids 31, 577 (1988).
• von Goeler, Stodiek, and Sauthoff (1974) S. von Goeler, W. Stodiek,  and N. Sauthoff, Phys. Rev. Lett. 33, 1201 (1974).
• Huysmans (2005) G. T. A. Huysmans, Plasma Phys. Contr. F. 47, B165 (2005).
• Candy and Waltz (2003) J. Candy and R. E. Waltz, J. Comput. Phys. 186, 545 (2003).
• Candy (2005) J. Candy, Phys. Plasmas 12, 072307 (2005).
• Kotschenreuther, Rewoldt, and Tang (1995) M. Kotschenreuther, G. Rewoldt,  and W. Tang, Comput. Phys. Commun. 88, 128 (1995).
• Pueschel, Kammerer, and Jenko (2008) M. J. Pueschel, M. Kammerer,  and F. Jenko, Phys. Plasmas 15, 102310 (2008).
• Parker, Chen, and Kim (2000) S. E. Parker, Y. Chen,  and C. C. Kim, Comput. Phys. Commun. 127, 59 (2000).
• Chen and Parker (2001) Y. Chen and S. Parker, Phys. Plasmas 8, 441 (2001).
• Chen and Parker (2003) Y. Chen and S. E. Parker, J. Comput. Phys. 189, 463 (2003).
• Cummings (1995) J. C. Cummings, Gyrokinetic simulation of finite-beta and self-generated sheared-flow effects on pressure-gradient-driven instabilities, Phd thesis, Princeton University (1995).
• Mishchenko, Hatzky, and Könies (2004) A. Mishchenko, R. Hatzky,  and A. Könies, Phys. Plasmas 11, 5480 (2004).
• Hatzky, Könies, and Mishchenko (2007) R. Hatzky, A. Könies,  and A. Mishchenko, J. Comput. Phys. 225, 568 (2007).
• Mishchenko et al. (2014a) A. Mishchenko, M. Cole, R. Kleiber,  and A. Könies, Phys. Plasmas 21, 052113 (2014a).
• Mishchenko et al. (2014b) A. Mishchenko, A. Könies, R. Kleiber,  and M. Cole, Phys. Plasmas 21, 092110 (2014b).
• Lee et al. (2001) W. W. Lee, J. L. V. Lewandowski, T. S. Hahm,  and Z. Lin, Phys. Plasmas 8, 4435 (2001).
• Startsev et al. (2016) E. Startsev, W.-L. Lee, W. Wang,  and Z. Lu, in B. Am. Phys. Soc. (American Physical Society, 2016) p. BM9.00002.
• Chen and Chacón (2014) G. Chen and L. Chacón, Comput. Phys. Commun. 185, 2391 (2014).
• Chen and Chacón (2015) G. Chen and L. Chacón, Comput. Phys. Commun. 197, 73 (2015).
• Chacón and Chen (2016) L. Chacón and G. Chen, J. Comput. Phys. 316, 578 (2016).
• Dudson et al. (2009) B. Dudson, M. Umansky, X. Xu, P. Snyder,  and H. Wilson, Comput. Phys. Commun. 180, 1467 (2009).
• Liu et al. (2014) Z. X. Liu, X. Q. Xu, X. Gao, T. Y. Xia, I. Joseph, W. H. Meyer, S. C. Liu, G. S. Xu, L. M. Shao, S. Y. Ding, G. Q. Li,  and J. G. Li, Phys. Plasmas 21, 090705 (2014).
• Fu et al. (2006) G. Y. Fu, W. Park, H. R. Strauss, J. Breslau, J. Chen, S. Jardin,  and L. E. Sugiyama, Phys. Plasmas 13, 052517 (2006).
• Breslau, Park, and Jardin (2006) J. A. Breslau, W. Park,  and S. C. Jardin, Journal of Physics: Conference Series 46, 97 (2006).
• Ferraro and Jardin (2009) N. Ferraro and S. Jardin, J. Comput. Phys. 228, 7742 (2009).
• Ren et al. (2015) X. Ren, M. Chen, X. Chen, C. Domier, N. Ferraro, G. Kramer, N. L. Jr., C. Muscatello, R. Nazikian, L. Shi, B. Tobias,  and E. Valeo, Journal of Instrumentation 10, P10036 (2015).
• Czarny and Huysmans (2008) O. Czarny and G. Huysmans, J. Comput. Phys. 227, 7423 (2008).
• Huysmans et al. (2010) G. Huysmans, S. Pamela, M. Beurskens, M. Becoulet,  and E. van der Plas, in Proceedings of the 23rd IAEA Fusion Energy Conference, Daejeon, South Korea (2010) pp. 11–16.
• Ku et al. (2016b) S. Ku, R. Hager, C. Chang, J. Kwon,  and S. Parker, J. Comput. Phys. 315, 467 (2016b).
• Chen et al. (2015) Y. Chen, J. Chowdhury, S. E. Parker,  and W. Wan, Phys. Plasmas 22, 042111 (2015).
• Balay et al. (2016a) S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang,  and H. Zhang, “PETSc Web page,” http://www.mcs.anl.gov/petsc (2016a).
• Balay et al. (2016b) S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, H. Zhang,  and H. Zhang, “PETSc users manual,” Tech. Rep. ANL-95/11 - Revision 3.7 (Argonne National Laboratory, 2016).
• Balay et al. (1997) S. Balay, W. D. Gropp, L. C. McInnes,  and B. F. Smith, in Modern Software Tools in Scientific Computing, edited by E. Arge, A. M. Bruaset,  and H. P. Langtangen (Birkhäuser Press, 1997) pp. 163–202.
• Cai and Fu (2012) H. Cai and G. Fu, Phys. Plasmas 19, 072506 (2012).
• Guazzotto et al. (2004) L. Guazzotto, R. Betti, J. Manickam,  and S. Kaye, Phys. Plasmas (1994-present) 11, 604 (2004).