# Computational Methods in Plasma Physics 2019¶

These are notes, lecture slides and code for use in PPPL’s Graduate Summer School, Computational Methods in Plasma Physics mini-course. All material is copyrighted by Ammar Hakim and released under the Creative Commons CC BY License.

Contents

- Computational Methods in Plasma Physics 2019
- Background
- Github repo for lectures and code
- Lecture 1: Introduction and Specialized Ordinary Differential Equations Solvers
- Lecture 2: The Boris algorithm and FDTD and FV schemes for Maxwell equations
- Lecture 3: The cutting-edge of computational plasma physics research
- Lecture 4: Discontinuous Galerkin Schemes
- References

## Background¶

Vast majority of plasma physics is contained in the Vlasov-Maxwell equations that describes the evolution of a particle distribution \(f_s(t,\mathbf{x},\mathbf{v})\) function in 6D phase-space. The particles move in electromagnetic fields that come from two sources: external coils and electrodes and the fields generated by the motion of the particles themselves. This makes the problem highly nonlinear: the fields tell the particles how to move, and the particles motion itself modifies the fields.

where \(\mathbf{a}_s = q_s/m_s(\mathbf{E}+\mathbf{v}\times\mathbf{B})\) is the acceleration due to Lorentz force on the particle and the right-hand side describes the effect of collisions. Of course, the electromagnetic fields are determined from Maxwell equations

Suitable modifications are required to account for relativistic effects.

Having written the basic equations it may appear that the field of plasma physics is complete. However, this is far from the case. The Vlasov-Maxwell system is a formidable system of coupled, nonlinear equations and describes vast physics that spans an enormous range of temporal and spatial scales. Everything from electron oscillations to slow, resistive evolution of objects in near equilibrium is contained within them. Hence, over the decades many approximations to the Vlasov-Maxwell equations have been developed that are often more tractable, allowing one to obtain reasonable results in specific physical situations. Of course, the frontier in computational and theoretical plasma physics remains the complete kinetic understanding of plasma from the full Vlasov-Maxwell equations, and much of the recent research focus has been on solving them (or asymptotic approximations in strong magnetic fields) directly.

In this mini-course we will look at the key modern computational techniques to handle various problems in plasma physics. This field remains an active area of research, and in the recent years has undergone a renaissance of sorts. My goal will be give you a flavor of the methods and point out key literature in this field. In particular we will study:

- Near first-principles simulations methods in which the Vlasov-Maxwell equation is solved directly. We will study two methods: the finite-difference-time-domain particle-in-cell method (FDTD PIC) (Lectures 1 and 2)
- Solving fluid equations that are obtained from taking moments of the Vlasov-Maxwell equations and making a closure to approximate the moments not evolved by the fluid equations. Here we will study the modern approach based on the theory of hyperbolic PDEs and Riemann solvers. We will also look at the special requirements of solvers that are required to study fusion problems (Lecture 3).
- Directly discretizing the Vlasov-Maxwell equations as a PDE in 6D. This is an emerging area of active research and has applications to study of turbulence in fusion machines and also exploring fundamental plasma physics in phase-space. (Lecture 4)

Through these topics we will look at the general applied mathematics literature for solvers for ordinary differential equations (ODEs) and partial differential equations (PDEs), with emphasis of those techniques that are useful in plasma physics.

## Github repo for lectures and code¶

You can get all the tex files for the slides, source for this website and the code I ran in class by cloning the Github repo for this class.

To run the code you need to install the Gkeyll software. Please see the Gkeyll website for instructions on how to get the code and the postprocessing tools needed. In brief:

Install the Python (mini)conda package manager

Install the Gkeyll code using:

conda install -c gkyl gkeyll

Install the postprocessing tool as:

conda install -c gkyl postgkyl

## Lecture 1: Introduction and Specialized Ordinary Differential Equations Solvers¶

### Code¶

See files code/lec1/sho-fwd-euler.lua and code/lec1/sho-mid-point.lua for forward Euler and mid-point schemes for harmonic oscillator problem. You can run the code from the directory it is kept as:

```
gkyl sho-mid-point.lua
```

To plot the trajectory do:

```
pgkyl -f sho-mid-point_ptclData.bp traj -e 90 -a 0.0 --fix-aspect
```

You can get help for options for the traj command by doing:

```
pgkyl traj --help
```

To plot the exact trajectory and the computed trajectory do:

```
pgkyl -f sho-mid-point_ptclData.bp -f sho-mid-point_exactData.bp traj -e 90 -a 0.0 --fix-aspect
```

You can save the animate to an mp4 file by passing the –save option to the traj command. This requires that you have the ffmpeg package installed.

### Summary¶

Read above introductory summary for overall background. It is important to understand how to derive conservation and other properties of the Vlasov-Maxwell system. Much of modern computational plasma physics is focused on inventing schemes that preserve at least some of these properties of the continuous system. See any plasma text book or the following excerpt from [Juno2019] for proofs.

The concept of phase-space volume preserving and symplectic schemes can be more easily understood by looking at the example of a simple harmonic oscillator

where \(\omega\) is the oscillation frequency.

To fully understand the physics behind these concepts one needs to
understand the *Lagrangian and Hamiltonian formulation* of
mechanics. For example, see text book of Goldstein
or first volume of Landau and Lifshitz, Mechanics. An overview of
Hamiltonian mechanics using *noncanonical coordinates* as applied to
single particle motion is given in Section II of [CaryBrizard2009].

A good description of various ODE solvers and their properties is given in Chapter 2 of [DurranBook]. Also see for formulas of the Strong-Stability preserving RK methods and their stability regions.

Several ODE schemes have been designed to handle stiff sources and in particular, diffusion terms arising from discretization of diffusion equations. See [Abdulle2013] and also [Meyer2013] for description of these schemes. In particular, the scheme by Meyer at al is to be preferred to it superior stability properties.

The ODE solvers described above are low order, that is second or third order. Some recent work attempts to construct very high order schemes (10-15th order!) that essentially makes the issues of conservation and other numerical errors mostly moot. For example, see [ReinSpiegel] for a 15th order scheme for use in gravitational N-body simulations. Such very high-order schemes have not found use in plasma-physics yet, mainly as the Maxwell solvers used in PIC codes are mostly second-order anyway. However, it is possible that these very high-order methods are useful in orbit codes.

## Lecture 2: The Boris algorithm and FDTD and FV schemes for Maxwell equations¶

PDF of Lecture 2 slides. Solution to the problem of finding \(\mathbf{A}\) if \(\mathbf{A} = \mathbf{R} + \mathbf{A}\times\mathbf{B}\) is here.

### Code¶

See files code/lec2/lorentz-boris.lua

You can play with this file to do various static or time-dependent electromagnetic fields. For example, motion in a constant magnetic field, in a field with a gradient, and in a driven system. See field specification in this write up for both non-resonant and resonant drivers.

### Summary¶

Particle-in-cell methods are based on pushing macro-particles. These represent the motion of characteristics in phase-space, along which the distribution function is conserved. The macro-particle equations-of-motion are

The most widely used method to solve this system of ODEs is the *Boris
algorithm*. See this excerpt from Birdsall and Langdon
book for details on how to implement this efficiently.

The Boris algorithm is surprisingly good: it is a *second-order*,
*time-centered* method that *conserves phase-space volume*. However,
the error in phase-velocity (that is there is an error in time-period
of orbits) accumulates *linearly*, as we saw for the harmonic
oscillator. See [Qin2013] for proofs that the Boris algorithm is
*not* symplectic but conserves phase-space volume.

The relativistic Boris algorithm does not compute the correct \(\mathbf{E}\times\mathbf{B}\) velocity. This can be corrected for and still maintain the volume-preserving property and was done in [HigueraCary2017].

The Yee-cell preserves the underlying geometric structure of Maxwell equations, and ensures that the divergence relations are maintained in the case of vacuum fields. In a plasma, however, current deposition needs to be done carefully to ensure current continuity is satisfied. See [Esirkepov2001], for example.

For extension of standard FDTD method to complex geometries, see, for example [Nieter2009] and other references. Recent research has focused on developing finite-element based PIC codes (that maintain geometric structure of Maxwell equations), but these are usually very expensive to run and very complex to develop.

Sometimes finite-volume schemes are also used to solve Maxwell equations. These may have some advantages and disadvantages compared to standard FDTD schemes. For example, FV usually do not conserve energy and find it hard to satisfy divergence relations. For a comparison of FV and FDTD methods see this page.

A comprehensive review of structure preserving algorithms for use in plasma physics is provided by [Morrison2017]. It has numerous references to the literature and should be consulted to develop a detailed understanding of such schemes.

## Lecture 3: The cutting-edge of computational plasma physics research¶

I showed slides from Eric Shi’s thesis defense. Please see his thesis for full details on the gyrokinetic algorithm, validation with LAPD and also first results from simulating the NSTX scrap-off-layer plasma.

### The Big Picture¶

Originally, I had intended to spend time talking about fluid solvers in this lecture. However, based on the extensive questions in the last lecture I decided the time would be better spent in looking at a few big-picture issues. We will return to numerics in Lecture 4.

In particular:

- What are the cutting-edge research questions in computational plasma physics?
- What is the relationship between modern numerical methods and
experiments and observations? (That is, why care about this stuff in
the first place? Can simulations
*predict*rather than*postdict*?) - How to incorporate “real-world” effects into simulations? (For example, boundary conditions, atomic physics, etc)

One can look at computational physics in two ways: as an end in itself, and as a tool for applications. Both of these are important!

As an end in itself:

- The first sits between applied mathematics and theoretical physics. The goal is to design efficient numerical methods to solve equations from theoretical physics.
- The goal here is the numerical method itself: what are its properties? Does it faithfully represent the underlying physics? Does it run efficiently on modern computers? Research into modern numerical methods (including structure preserving methods) fall into this category.
- Usually, besides the fun of solving complex equations (and writing
code), the goal is to gain deeper understanding of underlying
physics.
**Some theoretical questions can only be answered with computer simulations.** - This is a perfectly legitimate research area even if no connection to experiments is made, but only satisfies the curiosity of the researchers and helps one gain a better understanding of the physics.

As a tool for applications:

- The second is to look at the computational physics as providing tools to understand/design experiments or observations.
- Note that a large number of routine calculations are needed to build
modern experiments (heat-transfer, structural analysis, basic fluid
mechanics, equilibrium and stability calculations, etc).
**Such routine calculations are no longer cutting edge research topics**.

At the intersection of cutting-edge computational physics and modern
plasma physics is a set of **Billion Dollar Questions**. (In general,
one should not put currency values to such things).

These **Billion Dollar Questions** need huge investments in
experimental and observational programs as well as the very latest in
computational physics research.

Space Physics Examples: Parker Solar Probe

- Paker Solar Probe. “The
primary science goals for the mission are to
**trace how energy and heat move through the solar corona**and to explore**what accelerates the solar wind**as well as**solar energetic particles**.” - The Probe will collect detailed measurements of electric and magnetic fields as well as detailed distribution functions of particles.
- The solar wind plasma is nearly collisionless. It is likely that a proper understanding of kinetic physics (at the level of the Vlasov-Maxwell equations) will be needed to fully understand the physical processes.
- Cutting-edge simulations will be critcial to this. Serious research into numerics of Vlasov-Maxwell needs to be done and very large simulations need to be run.

Many other missions are active and planned: BepiColombo to Mercury; Juno to Jupiter.

- Much of the deep understanding of plasma processes in solar system
planets (magnetospheres, ionosphere) can only be gained from
detailed modeling: global kinetic modelling is likely
impossible.
**How to incorporate some kinetic effects into fluid models?**

Fusion Physics Examples: Building a working thermonuclear fusion reactor.

- The Iter project aims to build the world’s largest tokamak, a “magnetic bottle” to contain super-hot plasma and heat it to ignition temperatures.

There are other major fusion efforts around the world:

- Beautiful stellarators (and Wiki article) that may have better properties than tokamaks and provide a faster route to fusion energy
- High-field based compact tokamaks; field-reversed configurations; spinning magnetic mirror machines; etc

There are major unsolved problems in the basic physics of fusion machines. Most of these can only be answered by large-scale computing and much of the numerical tools have not yet been fully developed.

The Scientific Discovery through Advanced Computing program in
fusion has large projects that address the very serious **Billion
Dollar Question**: will controlled fusion be eventually possible?

- The numerics research here is focussed on gyrokinetic and even full
kinetic understanding of fundamental turbulence and transport
processes in the tokamak.
**These equations are very difficult to solve!** - Disruptions are dangerous processes that can “kill” certain fusion machine: large-scale MHD simulations are needed. Significant new research is being done in new numerical methods and application of existing MHD codes to such problems.
- Runaway electrons (relativistic high-energy electron beams) can drill holes in fusion machines. See SCREAM project and special PPCF issue.
- Very serious!
**Will need huge kinetic calculations**. Also, the formulation of self-consistent coupling betwen the runaway electrons and MHD is not complete. See review by [Boozer2015].

These are only selection of problems I am directly familiar with. I hope it gives you a flavor and understanding why computational plasma physics is such a serious and important field!

## Lecture 4: Discontinuous Galerkin Schemes¶

### Summary¶

The discontinuous Galerkin method is a type of finite-element scheme
in which the solution across cell boundaries can be discontinuous. The
method was originally invented for elliptic equations by Nitsche in
1971 (paper is in German). However, the key paper on application to a
hyperbolic PDE was written by Reed and Hill in 1973. The latter paper
has more than 2100 citations. A later key paper on extension of the
method to nonlinear systems of equations was by Cockburn and Shu
(JCP **41**, 199-224 1998).

The key ingredient of DG is Galerkin minimization on a
finite-dimensional *discontinuous* space. This space is usually
piecewise continuous polynomials, but could be other types of
functions too.

I introduced the idea of weak-equality, which is central to a proper understanding of the DG scheme. This weak-equality concept, and idea of recovery built from it, can be used to construct schemes for both advection and diffusion equation that is higher-order accurate than standard DG in smooth regions. Some form of limited recovery is needed when the solution is not smooth.

One can solve the full Vlasov-Maxwell equations using DG. See
[Juno2019] for example. The scheme described there conserves energy
exactly, but does not conserve momentum. This is a typical feature of
solvers for Vlasov equations: one can *either* conserve momentum or
energy but not both. The construction of a momentum *and* energy
conserving scheme remains an ongoing research problem.

## References¶

[Juno2019] | (1, 2) Juno, J., Hakim, A., TenBarge, J., Shi, E.,
Dorland, W. (2018). “Discontinuous Galerkin algorithms for fully
kinetic plasmas”, Journal of Computational Physics, 353,
110–147. http://doi.org/10.1016/j.jcp.2017.10.009 |

[CaryBrizard2009] | Cary, J. R., &
Brizard, A. J. (2009). “Hamiltonian theory of guiding-center
motion”. Reviews of Modern Physics, 81 (2),
693–738. http://doi.org/10.1103/RevModPhys.81.693 |

[DurranBook] | Dale E. Durran, “Numerical Methods for Fluid Dynamics”, Springer. Second Edition. |

[Abdulle2013] | Abdulle, A., & Vilmart, G. (2013). “PIROCK: A
swiss-knife partitioned implicit–explicit orthogonal Runge–Kutta
Chebyshev integrator for stiff diffusion–advection–reaction
problems with or without noise”. Journal of Computational Physics,
242 (C), 869–888. http://doi.org/10.1016/j.jcp.2013.02.009 |

[Meyer2013] | Meyer, C. D., Balsara, D. S., & Aslam, T. D. (2014). “A
stabilized Runge–Kutta–Legendre method for explicit
super-time-stepping of parabolic and mixed equations”. Journal of
Computational Physics, 257 (PA),
594–626. http://doi.org/10.1016/j.jcp.2013.08.021 |

[ReinSpiegel] | Rein, H., & Spiegel, D. S. (2014). ias15: a fast, adaptive, high-order integrator for gravitational dynamics, accurate to machine precision over a billion orbits. Monthly Notices of the Royal Astronomical Society, 446(2), 1424–1437. http://doi.org/10.1093/mnras/stu2164 |

[Qin2013] | Qin, H., Zhang, S., Xiao, J., Liu, J., Sun, Y., &
Tang, W. M. (2013). “Why is Boris algorithm so good?” Physics of
Plasmas, 20 (8), 084503–5. http://doi.org/10.1063/1.4818428 |

[HigueraCary2017] | Higuera, A. V., &
Cary, J. R. (2017). “Structure-preserving second-order integration
of relativistic charged particle trajectories in electromagnetic
fields”. Physics of Plasmas, 24 (5),
052104–7. http://doi.org/10.1063/1.4979989 |

[Esirkepov2001] | Esirkepov, T. Z. (2001). “Exact charge conservation
scheme for Particle-in-Cell simulation with an arbitrary
form-factor”, Computer Physics Communications, 135, 144–153. |

[Nieter2009] | Nieter, C., Cary, J. R., Werner, G. R., Smithe, D. N.,
& Stoltz, P. H. (2009). “Application of Dey–Mittra conformal
boundary algorithm to 3D electromagnetic modeling”. Journal of
Computational Physics, 228 (21),
7902–7916. http://doi.org/10.1016/j.jcp.2009.07.025 |

[Morrison2017] | Morrison, P. J. (2017). Structure and
structure-preserving algorithms for plasma physics. Physics of
Plasmas, 24 (5), 055502–21. http://doi.org/10.1063/1.4982054 |

[Boozer2015] | Boozer, A. H. (2015). “Theory of runaway electrons in
ITER: Equations, important parameters, and implications for
mitigation”. Physics of Plasmas, 22 (3),
032504–18. http://doi.org/10.1063/1.4913582 |