2 Physical model

As already mentioned in the introduction we will consider a kinetic model of a quasi-neutral edge plasma. Thus, we start with the inhomogeneous Vlasov equation (for the electrons)

where the electron density $f$ is a function of position $\boldsymbol{r}$ and velocity $\boldsymbol{v}$ , and the electron mass is denoted by $m_{e}$ . The collision term $C(f)$ describes the particle–particle interactions.

In our model of the edge plasma we transform the problem into a coordinate system that is aligned to the magnetic field, which consequently is parallel to the $z$ -axis. This is the same coordinate system that is used in Chodura (Reference Chodura1992). By assuming a homogeneous distribution in the directions perpendicular to the magnetic field we eliminate all but one space coordinate. As a matter of fact, the model itself does not use the magnetic field further. Thus the electric field $\boldsymbol{E}$ is the only field introduced into the equations and consequently the electrons experience the force of $\boldsymbol{F}=(0,0,eE_{z})$ , where $e$ is the charge and $E_{z}$ is the $z$ -component of the electric field. We assume that all other components of the electric field are negligible. This assumption is essential for the coordinate transformation we employ later on. In the present work we further assume that the electron and ion temperatures are of the same order of magnitude. Therefore the velocity of the ions is much slower than the velocity of the electrons. This implies that we can neglect the ion velocity in the electron–ion collision term. We also assume that $E_{z}$ is constant (in the present model this is a good approximation; see for example Stangeby Reference Stangeby2000; Tskhakaya, Groth & JET EFDA Contributors Reference Tskhakaya and Groth2015).

We further assume that the electrons are perfectly magnetized and the corresponding distribution function depends on the parallel ( $v_{z}$ ) and total ( $v$ ) velocities. For convenience, below we introduce the cosine $\unicode[STIX]{x1D707}=\cos \unicode[STIX]{x1D703}=v_{z}/v$ of the pitch angle $\unicode[STIX]{x1D703}$ . In this coordinate system the Vlasov equation takes the form

Although, for the plasma edge, the collision operator in general includes Coulomb interactions (see, for example, Karney Reference Karney1986) as well as charged–neutral (for example, we refer to Janev, Reiter & Samm (Reference Janev, Reiter and Samm2004)) and inelastic charged–charged particle collisions (see, for example, Tskhakaya et al. (Reference Tskhakaya and Groth2015), and the references therein), in the present work we restrict ourselves to Coulomb electron–ion collisions. Due to the fact that the ions are much slower than the electrons we use the pitch-angle collision operator, which is a simplification of the operator derived in Rosenbluth, MacDonald & Judd (Reference Rosenbluth, MacDonald and Judd1957):

where $\unicode[STIX]{x1D6FE}_{ei}$ is a constant (in our model) comprising all physically relevant parameters

This formula introduces the electron density $n_{e}$ , the Landau logarithm $L$ , the vacuum permittivity $\unicode[STIX]{x1D716}_{0}$ and the thermal speed $v_{T}$ . The specific values are listed in table 1.

Table 1.

All used parameters are shown. They represent typical values for an edge plasma in a fusion device. The employed coefficients $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D706}$ are also listed.

For the numerical simulation we render (2.2) and (2.3) dimensionless. Thus, we define the velocity normalized by the thermal speed

and measure length as

where $Z$ is the length of the spatial domain. These two quantities imply the following normalization of time

Note that no normalization is required for $\unicode[STIX]{x1D707}$ as this is already a dimensionless variable.

Finally, the dimensionless equation used in the numerical simulation is as follows

Since we are only concerned with steady state solutions, we drop the time derivative in the present paper and obtain

Combining all the relevant parameters we finally obtain

where we have defined

and

Typical values for the dimensionless parameters in edge plasmas can be found in table 1.

From a physical point of view (2.10) has three contributions. The first term on the left-hand side describes an advection along the magnetic field lines with speed $\unicode[STIX]{x1D707}\,\unicode[STIX]{x1D709}=v_{z}/v_{T}$ . The acceleration due to the static electric field is modelled by the second term on the left-hand side. Notice that this results in an advection for $\unicode[STIX]{x1D707}$ as the electric field along the $z$ -axis increases/decreases the velocity along $z$ and thus increases/decreases $\unicode[STIX]{x1D707}$ . Finally, the right-hand side describes the diffusive effects due to the slow ions pitch-angle scattering the electrons.

All the terms in (2.10) scale differently with respect to $\unicode[STIX]{x1D709}$ . For large $\unicode[STIX]{x1D709}$ the $\unicode[STIX]{x1D701}$ -advection is prevalent. From a physical point of view this means that the electrons are sufficiently fast such that they do not experience many interactions with the ions. Thus diffusive effects are weak. In contrast, for small $\unicode[STIX]{x1D709}$ , the slow electrons experience a strong diffusion. Note that the used collision term is not valid for very small $\unicode[STIX]{x1D709}$ , since the electrons are assumed to be much faster than the ions.

To set the boundary condition for the above problem we have to keep in mind that we want to describe a realistic domain for edge plasmas, i.e. plasma as it can be observed e.g. in the SOL. Due to the assumed symmetry in the angle-like coordinate $\unicode[STIX]{x1D707}$ we infer that any electron flux crossing the $\unicode[STIX]{x1D707}$ boundary must be balanced by a flux from the other side of the boundary. Thus, we will use homogeneous Neumann boundary conditions in this case. For the $\unicode[STIX]{x1D701}$ -boundary we have to make a distinction between the left and the right boundary. On the one hand, at the boundary at $\unicode[STIX]{x1D701}=0$ , the inflow is determined by the relevant plasma dynamics in the interior of the device. The actual value is highly problem dependent, but we generally use Dirichlet boundary conditions there. On the other hand at the right boundary at $\unicode[STIX]{x1D701}=1$ we have to model the wall of the device. As a first approximation we use homogeneous Dirichlet boundary conditions there. In other words we consider the case where the magnetic field is normal to the wall and the electrons are magnetized in the entire simulation domain. However, more involved boundary conditions that reflect a certain amount of plasma are certainly a viable (and probably more realistic) alternative that we will explore in the future. The boundary configurations used in the following examples are illustrated in figure 1.

Figure 1.

In this illustration the two-dimensional grid in $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D707}$ is shown with the imposed boundary conditions. Boundaries with homogeneous Neumann boundary conditions are coloured blue, whereas those with Dirichlet boundary conditions are coloured red. Uncoloured boundaries are outflow boundaries. Next to the inflow boundaries two possible boundary data are sketched. These isotropic and electron beam inflow conditions are tagged with

and 

, respectively.

3 Numerical methods

In this section we will describe the numerical methods used to obtain an efficient approximation to (2.10). To that end we will first describe the discretization and then discuss the iterative methods that are used to solve the resulting linear system. The latter is particularly important in order to obtain an efficient numerical method.

Let us note that $f$ depends on $\unicode[STIX]{x1D701}$ , $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D709}$ . However in the present case $\unicode[STIX]{x1D709}$ is a parameter only and we thus have to solve a two-dimensional problem in $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D707}$ repeatedly (once for each value of $\unicode[STIX]{x1D709}$ ). It should be kept in mind, however, that depending on the value of $\unicode[STIX]{x1D709}$ the relative strength of diffusion and advection can differ significantly.

First, the diffusive term is discretized by the following finite difference approximation

In the last equation we denoted the equidistant $\unicode[STIX]{x1D707}$ -step size by $h_{\unicode[STIX]{x1D707}}$ and $f_{i}=f(\unicode[STIX]{x1D701},\unicode[STIX]{x1D707}_{i})$ , where $\unicode[STIX]{x1D707}_{i}=ih_{\unicode[STIX]{x1D707}}$ for $i\in \{-n_{\unicode[STIX]{x1D707}},\ldots ,n_{\unicode[STIX]{x1D707}}\}$ with $n_{\unicode[STIX]{x1D707}}h_{\unicode[STIX]{x1D707}}=1$ .

In principle, any second-order accurate stencil can be used here. Note, however, that the stencil above results in a symmetric matrix. As has been pointed out in the previous section, Neumann boundary conditions are used for both the boundary at $\unicode[STIX]{x1D707}=1$ and the boundary at $\unicode[STIX]{x1D707}=-1$ . The additional $\unicode[STIX]{x1D707}$ -advection term, which is the second term on the left-hand side of (2.10), is discretized with centred finite differences. Note, however, that this term degenerates and thus no boundary conditions are required.

In contrast, for the $\unicode[STIX]{x1D701}$ -advection we use an upwind scheme (i.e. a different stencil is employed for positive and negative $\unicode[STIX]{x1D707}$ , respectively). This facilitates the implementation of the Dirichlet boundary conditions and guarantees that we have perfect outflow conditions in the regions where this is important (we refer again to figure 1).

The resulting linear system can be solved with an arbitrary numerical method. However, our goal is to develop an efficient approach. Although we only consider the $1+2$ dimensional case in the present paper, it is an important consideration that such schemes can be generalized to the higher-dimensional setting. In this context it is of particular interest to discuss numerical methods that scale well with respect to memory consumption. Consequently, we will employ so-called iterative methods, which can be implemented matrix free. That is, only the application of the matrix to a vector has to be computed (and there is no need to assemble and store the entire matrix in memory). The disadvantage of these methods, however, is that they usually require a large number of iterations which results in a slow algorithm. In many applications this problem can be overcome by employing a suitable preconditioner.

In the described linear problem we have a matrix $\unicode[STIX]{x1D63C}$ and a right-hand side $\boldsymbol{b}$ . In the present case $\boldsymbol{b}$ has very few non-zero elements, which are due to the boundary terms. Our goal is then to solve

by some iterative method (see, for example, Saad Reference Saad2003).

A preconditioner $\unicode[STIX]{x1D648}$ is a matrix that can be easily inverted and constitutes an approximation to $\unicode[STIX]{x1D63C}$ . For a right preconditioner we set $\widetilde{\unicode[STIX]{x1D63C}}=\unicode[STIX]{x1D63C}\unicode[STIX]{x1D648}^{-1}$ , $\boldsymbol{y}=\unicode[STIX]{x1D648}\boldsymbol{x}$ and consider the linear system

which has to be solved for $\boldsymbol{y}$ . The solution of this preconditioned system hopefully requires less iterations as ideally $\unicode[STIX]{x1D63C}\unicode[STIX]{x1D648}^{-1}$ is close to the identity. The unknown $\boldsymbol{x}$ can then be determined easily by solving $\unicode[STIX]{x1D648}\boldsymbol{x}=\boldsymbol{y}$ .

Similar to this approach a left preconditioner takes the form

This, however, has the disadvantage that the residual is modified which causes some complications in deciding when to stop the iteration. Consequently, we will only use right preconditioners in the present work.

Let us now turn our attention to constructing an efficient preconditioner. For elliptic and parabolic problems preconditioners based on multigrid methods (Briggs, Henson & McCormick Reference Briggs, Henson and McCormick2000) or domain decomposition (Smith, Bjorstad & Gropp Reference Smith, Bjorstad and Gropp2004) can be used. While some of these methods have been generalized to hyperbolic problems (in particular, in the context of fluid models; see, for example, Adams, Samtaney & Brandt Reference Adams, Samtaney and Brandt2010; Shadid et al. Reference Shadid, Pawlowski, Banks, Chacónc, Lin and Tuminaro2010) significant challenges remain. Consequently, in recent years so-called physics-based preconditioners have become popular (see, for example, Chacón & Knoll Reference Chacón and Knoll2003; Chacón Reference Chacón2008; Reynolds, Samtaney & Woodward Reference Reynolds, Samtaney and Woodward2010; Reynolds, Samtaney & Tiedeman Reference Reynolds, Samtaney and Tiedeman2012). These preconditioners exploit the particular structure of the equation under consideration and are often successful where more generic algorithms are inefficient or fail altogether.

In our problem, the matrix $\unicode[STIX]{x1D63C}$ can be partitioned into terms that only contain derivatives in the $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D701}$ directions, respectively. Thus, we can write our matrix as $\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D701}}+\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ . This partitioning is the basis of all the preconditioners developed in this paper. The linear system corresponding to the matrices $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D701}}$ and $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ can be solved efficiently (using, for example, the Thomas algorithm) as these are only one-dimensional problems. Put differently, we consider preconditioners based on dimension splitting.

As a first step, one might argue that the diffusion found in $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ is what limits the convergence speed. Thus, we define $\unicode[STIX]{x1D648}=\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ . In the following we refer to this as the $\unicode[STIX]{x1D707}$ -preconditioner. While this approach results in some improvement compared to the iterative schemes without any preconditioner the scaling is still far from optimal.

Therefore, we will consider a second preconditioner based on an ADI approach (such methods have a long tradition in the time integration of partial differential equations; see, for example, Yanenko (Reference Yanenko1971), Lindemuth & Killeen (Reference Lindemuth and Killeen1973), Schnack & Killeen (Reference Schnack and Killeen1980) or for more recent work Hujeirat (Reference Hujeirat1998), Sgura, Bozzini & Lacitignola (Reference Sgura, Bozzini and Lacitignola2012), Reynolds et al. (Reference Reynolds, Samtaney and Tiedeman2012), Müller & Scheichl (Reference Müller and Scheichl2014)). We define

where $s$ is a tuneable parameter. Note that this is equivalent to a Lie–Trotter splitting step in a time integration scheme (with step size $s$ ). Since performing a numerical time integration yields a successively better approximation to the steady state solution, we can use this splitting scheme as a preconditioner. This approach combines the advantages of iterative methods (which force a relaxation to the equilibrium of interest) with a time marching scheme (which allows us to construct an efficient preconditioner).

At first sight we might be tempted to choose a large $s$ (since an approximation at a later time yields a better approximation to the steady state solution). However, due to the splitting error committed there is a trade-off to be made here. In the numerical simulations conducted we found that values of the order of $0.1$ give the fastest convergence.

Both of these preconditioners can be implemented without storing the corresponding matrix and suitable $O(n)$ algorithms are available for each step. Thus, the overhead in computational cost is approximately a factor of three compared to the unpreconditioned case. As we will see in § 4 the ADI preconditioner reduces the number of iterations by two orders of magnitude. Consequently the algorithm based on the ADI preconditioner is significantly less computationally demanding.

Both of these preconditioners are generic enough that they can be combined with any iterative numerical method. In particular, we will consider Richardson iteration and the generalized minimal residual method (GMRES) (see, for example, Saad Reference Saad2003). The latter is part of the class of Krylov subspace methods and usually converges faster than relaxation schemes. This performance comes at the cost that we have to store multiple vectors of previous iterations in memory. In contrast, Richardson iteration can be implemented using a single input and a single output array and is thus extremely memory efficient.

To conclude this section, let us discuss an approach to find the steady state solution that is commonly used in the literature. Namely, to integrate the time-dependent equations of motion (2.10) sufficiently long in order to obtain an approximation to the steady state. However, in the present case this approach suffers from severe difficulties. This is due to the fact that the present model includes an advection in the $\unicode[STIX]{x1D707}$ direction the speed of which can be arbitrarily small (in the numerical discretization the smallest speed is proportional to the grid size $h$ ). Since we have to resolve the corresponding time scale, the final time has to scale as $1/h$ (which implies that the final time of the integration is large). On the other hand, performing splitting is attractive as it allows us to obtain an efficient time integrator without a Courant-Friedrichs-Lewy (CFL) constraint. Moreover, such methods have been studied extensively (see, for example, Cheng & Knorr Reference Cheng and Knorr1976; Klimas & Farrell Reference Klimas and Farrell1994; Sauter et al. Reference Sauter, Harvey and Hinton1994; Filbet et al. Reference Filbet, Sonnendrücker and Bertrand2001; Arber & Vann Reference Arber and Vann2002; Crouseilles et al. Reference Crouseilles, Mehrenberger and Sonnendrücker2010; Qiu & Christlieb Reference Qiu and Christlieb2010; Rossmanith & Seal Reference Rossmanith and Seal2011; Qiu & Shu Reference Qiu and Shu2011; Landreman & Ernst Reference Landreman and Ernst2012; Einkemmer & Ostermann Reference Einkemmer and Ostermann2014a ,Reference Einkemmer and Ostermann b ; Christlieb Reference Christlieb2015; Crouseilles et al. Reference Crouseilles, Einkemmer and Faou2015; Casas et al. Reference Casas, Crouseilles, Faou and Mehrenberger2017; Crouseilles et al. Reference Crouseilles, Einkemmer and Faou2016; Einkemmer Reference Einkemmer2016). However, the error made by the time integrator then dictates a relatively small time step (due to accuracy requirements) which implies that we have to take a large number of steps in order to integrate to the final time. Consequently, such methods are not computationally efficient for the problems considered here. Note that, in particular, the GMRES approach gives us more freedom in this regard as, in a manner of speaking, we can take different time steps for different parts of the problem.

4 Numerical results

To compare the different numerical methods we compute a solution to (2.10) for different grid sizes. We initialize the solution in the interior of the domain with zero and consider in this paper the following inflow boundary conditions:

1. (i) Isotropic particle distribution. The corresponding Dirichlet boundary conditions take the constant value $1$ as is illustrated in figure 1.

2. (ii) Electron beam. In this case we model an electron beam that enters the edge plasma but is not aligned to the magnetic field lines. For the corresponding Dirichlet boundary conditions a Gaussian beam profile given by

   (4.1) $$\begin{eqnarray}\displaystyle f(\unicode[STIX]{x1D707},\unicode[STIX]{x1D701})=\left\{\begin{array}{@{}ll@{}}\text{e}^{-100\cdot (\unicode[STIX]{x1D707}-0.5)^{2}},\quad & \unicode[STIX]{x1D707}\geqslant 0,\unicode[STIX]{x1D701}=0,\\ 0,\quad & \unicode[STIX]{x1D707}<0,\unicode[STIX]{x1D701}=1\end{array}\right. & & \displaystyle\end{eqnarray}$$
   is prescribed.

Numerical results for both of these boundary data are shown in figure 2.

Figure 2.

On (a) we see the numerically computed steady state when using isotropic Dirichlet boundary conditions for the inflow at $\unicode[STIX]{x1D701}=0,1$ . The results on (b) are obtained by using an electron beam instead. To create these plots we used $\unicode[STIX]{x1D705}=0.385$ , $\unicode[STIX]{x1D706}=0.0385$ , $\unicode[STIX]{x1D709}=0.5$ .

All algorithms stop when a certain precision is reached (in all the numerical simulations considered here we require that the residual is below $10^{-5}$ ). We only show results with an equal number of grid points in both the $\unicode[STIX]{x1D707}$ and $\unicode[STIX]{x1D701}$ direction. Thus, we can compare the number of iterations each method requires to achieve a given precision (averaged over a range of values for $\unicode[STIX]{x1D709}$ ). All errors are measured in the $L^{2}$ norm. The corresponding results are plotted in figure 3.

Figure 3.

The number of iterations each numerical method requires to reach a precision below $10^{-5}$ . These numbers are obtained by averaging the results of the parameters $\unicode[STIX]{x1D705}=0.385$ , $\unicode[STIX]{x1D706}=0.0385$ and $\unicode[STIX]{x1D709}=\{0.25,0.5,0.75,1.0,1.25,1.5,2.0,3.0,4.0,5.0,6.0\}$ . The parameter in the ADI scheme was set to $s=0.1$ . The ▫ and ○ respectively mark the results obtained with isotropic and beam initial values, as described in the text.

Figure 4.

The number of iterations each numerical method requires to reach a precision below $10^{-5}$ . These numbers are obtained by averaging the results of the parameters $\unicode[STIX]{x1D705}=0.385$ , $\unicode[STIX]{x1D706}=0.0385$ and $N=\{32,64,96,128,160,192,224,256,288,320,352\}$ . The parameter in the ADI scheme was set to $s=0.1$ . The ▫ and ○ respectively mark the results obtained with isotropic and beam initial values, as described in the text.

We can observe from figure 3 that the number of iterations required differs slightly with respect to the boundary condition used. We conclude that all preconditioners proposed in this paper improve the rate of convergence significantly. In addition, we find the following.

1. (i) Only inverting $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ is the cheapest preconditioner but also the least effective one. While this preconditioner might be competitive for small grids, as we increase the number of grids points it scales as $N^{2.2}$ , where $N$ is the number of grid points in one direction.

2. (ii) The ADI preconditioner works very well and scales almost ideally as we increase the number of grid points. As a result this preconditioner is able to reduce the number of iterations required by two orders of magnitude (compared to the unpreconditioned GMRES implementation). This is independent of the iterative method used (i.e. both GMRES as well as Richardson iteration show this behaviour). However, as expected, the GMRES implementation significantly outperforms the implementation based on Richardson iteration (by approximately a factor of $3$ ) (see Saad (Reference Saad2003) for instance). This comes at the cost of increased memory consumption (in the simulation we employ a restarted GMRES scheme with Krylov dimension $20$ ).

To conclude this section we investigate the speed of convergence as a function of $\unicode[STIX]{x1D709}$ . In figure 4 we see that, in general, the highest computational effort is necessary for small $\unicode[STIX]{x1D709}$ . Although it is for these $\unicode[STIX]{x1D709}$ that the GMRES scheme benefits the most from employing preconditioning. In addition, we observe that the relative difference in performance between the $\unicode[STIX]{x1D707}$ -preconditioner and the ADI preconditioner decreases for small $\unicode[STIX]{x1D709}$ . This is expected because in this limit the operator $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D707}}$ dominates $\unicode[STIX]{x1D63C}_{\unicode[STIX]{x1D701}}$ .