# Hermes: Global plasma edge fluid turbulence simulations

## Abstract

The transport of heat and particles in the relatively collisional edge regions of magnetically confined plasmas is a scientifically challenging and technologically important problem. Understanding and predicting this transport requires the self-consistent evolution of plasma fluctuations, global profiles and flows, but the numerical tools capable of doing this in realistic (diverted) geometry are only now being developed.

Here a 5-field reduced 2-fluid plasma model for the study of instabilities and turbulence in magnetised plasmas is presented, built on the BOUT++ framework. This cold ion model allows the evolution of global profiles, electric fields and flows on transport timescales, with flux-driven cross-field transport determined self-consistently by electromagnetic turbulence. Developments in the model formulation and numerical implementation are described, and simulations are performed in poloidally limited and diverted tokamak configurations.

###### pacs:

52.25.Xz, 52.65.Kj, 52.55.Fa^{1}

## 1 Introduction

The edge of magnetically confined plasmas, such as tokamaks, is where hot confined plasma encounters neutral gas, material surfaces, and the associated impurities. The transport of heat and particles in this region determines the heat loads and erosion rates of plasma facing components (PFCs). Predicting this transport, and exploring means of reducing heat fluxes to PFCs, has played an important role in designing the ITER divertor [1], and will be critical to the design of a future DEMO device [2, 3]. One of the uncertainties in making these predictions is the transport across the magnetic field, which is not well described by diffusion [4, 5], and is thought to be turbulent [6]. Since the fluctuations can be of similar spatial scales and magnitude to average profiles, significant effort has been devoted to developing and testing 3D flux-driven fluid turbulence simulation codes, including GBS [7, 8], TOKAM-3D [9] and TOKAM3X [10].

In this paper we present Hermes, a new tool for the simulation of collisional plasmas, in particular the edge and divertor regions of tokamaks. The ultimate aim of this work is to construct a model capable of simulating self-consistently the turbulence, plasma profiles and flows in the edge of magnetically confined plasmas over transport timescales. This is an ambitious undertaking, as it requires the combination of neutral gas and atomic physics, plasma-wall interaction, turbulent transport, and neoclassical effects [11]. As a first step towards this, we present in section 2 a drift-reduced model [12, 13, 14, 15] which has been constructed based on the derivation of Simakov and Catto [16], and similar to that derived in [17] but with several modifications to make it more suitable for numerical solution in global geometry. This model has good conservation properties, described in section 2.1, and has been implemented using the BOUT++ framework [18, 19] using conservative numerical methods described in section 3.

Significant advances have been made in improving numerical accuracy and stability. In particular, PETSc [20, 21, 22] has been used to implement a new solver for the axisymmetric potential , described in section 5.3. This allows BOUT++ simulations to self-consistently evolve the axisymmetric electric field in X-point geometry for the first time. Results are presented showing that important features of the global equilibrium can be recovered, including the Pfirsch-Schlüter current and Geodesic Acoustic Mode (GAM) oscillations [23].

## 2 Hermes model equations

The plasma model is electromagnetic, and evolves electron density , electron pressure where is the electron temperature, parallel ion momentum , plasma vorticity , and electromagnetic potential . There is no separation between background and fluctuations in this model, since we wish to study regions such as the tokamak edge where these are of similar magnitude. The model discussed here (Hermes-1) assumes cold ions; a hot-ion extension of this model (Hermes-2) is currently under development. The evolution equations normalised to the ion gyro-frequency and sound gyro-radius are:

(1) | |||||

(2) | |||||

(3) | |||||

(4) | |||||

(5) | |||||

with cross-field EB and magnetic drifts given by:

(6a) | |||||

(6b) |

Here we use the notation , , and . The electron beta appearing in equation 5 is . The parallel current is used to calculate the parallel electron velocity using . The parallel electron thermal conduction coefficient is the Braginskii value , where is the electron thermal speed and is the electron collision rate, with optional flux limiters as used in SOLPS [25]. The resistivity is given by . Anomalous diffusion can be represented in axisymmetric simulations (section 5.2) by particle and thermal diffusivities and .

The form of the diamagnetic current, in terms of a magnetic drift , is used here because it is more easily implemented in terms of fluxes through cell faces than the standard form, and hence suitable for the conservative numerical schemes described in section 3: The flow speed depends only on the local temperature, rather than on pressure gradients, and this has been found to improve numerical stability. This approach has also been used in the TOKAMX code [10].

The vorticity is related to the electrostatic potential by

(6g) |

where the Boussinesq approximation is used, replacing with a constant (with no spatial or temporal variation) in the vorticity equation. The capability to perform simulations without making this approximation has been used in BOUT++ for other models [26, 27] but is not used here because it cannot yet be used in diverted X-point geometry for the axisymmetric electric field (section 5.3). The calculation of electrostatic potential from the vorticity is a crucial component in all drift-reduced plasma models, including Hermes. Developments to allow equation 6g to be solved correctly in X-point geometry are described in section 5.3.

### 2.1 Conservation properties

The equations 1-5 are formulated in divergence form, and the operators are cast in terms of fluxes between cells to ensure conservation of the quantity being advected (see section 3). This is potentially important for advection of plasma density, since experience with transport codes has shown that results are sensitive conservation of particle number in high recycling regimes [1].

In order to study conservation of energy, the procedure described in [17, 28] is followed: Multiply the vorticity equation (3) by , and Ohm’s law (5) by . Rearranging, a set of divergence terms and a set of transfer channels are obtained. When integrated over the spatial domain the divergence terms become fluxes through the boundary, which can be made to vanish by appropriate choice of boundary conditions. The rate of change of each form of energy, in the absence of boundary fluxes or sources and neglecting cross-field diffusion and dissipation terms, is given by:

(6h) |

(6i) |

(6j) | |||||

(6k) | |||||

which correspond to the ion energy (equation 6h), ion parallel kinetic energy (equation 6i), electromagnetic field energy (equation 6j), and electron thermal energy (equation 6k). Each term on the right hand side of equations 6h-6k has a corresponding term in another equation with which it balances, representing a transfer of energy from one form to another. This model therefore has a conserved energy of the form

(6l) | |||||

Unfortunately this conservation is not complete in the implementation of this model, because the ion parallel momentum (equation 4) does not contain the ion polarisation drift. The first term on the right of equation 6i vanishes if the total ion velocity is used in the advection of ion parallel momentum equation, but in equation 4 the ion velocity is approximated by the sum of and parallel flow only (ion diamagnetic velocity being zero for cold ions). The ion polarisation drift is a higher order than the and diamagnetic drifts considered here, and including the ion polarisation drift in the ion momentum equation would introduce a time derivative into the right hand side, which could not then be solved by the Method of Lines approach adopted by BOUT++. This is a challenge of drift-reduced models, though a possible solution is suggested in [17] and will be investigated as part of future work.

### 2.2 Boundary conditions

At the radial boundaries simple Dirichlet or Neumann boundary conditions are used, chosen so that the flux of particles and energy through the boundary goes to zero. A difficulty arises at the core boundary, since the magnetic drift is approximately vertical, and so into the boundary at the top of the device and out of the boundary at the bottom (or vice-versa, depending on the sign of the drift). If this is artificially prevented from convecting particles and thermal energy through the core boundary then a narrow unphysical boundary layer forms. As a result here there can be a net flux of energy through the core boundary. In most cases this source is expected to be small compared to the external input power, but pathological cases could exist. A possible future improvement would be to ensure that the total energy and particle flux integrated over the core boundary vanishes.

In the direction parallel to the magnetic field sheath boundary conditions are needed. The correct boundary conditions to apply at the entrance to the sheath in magnetically confined plasmas is the subject of a long and ongoing debate [29, 30, 31, 32, 33]. Here we adopt relatively simple boundary conditions, leaving the choice of more complex boundary condition to future work. The parallel ion velocity is sonic:

(6m) |

where is the sound speed. In practice the inequality means that if the flow in front of the sheath becomes supersonic then the boundary condition becomes zero-gradient (Neumann condition). Allowing for this possibility is important, since several mechanisms can produce a supersonic transition, which have been studied analytically [34] and in simulations [35].

The parallel current at the sheath is given by:

(6n) |

Note that the ion speed into the sheath is used rather than just , to handle the supersonic regime correctly. The electron flow into the sheath saturates, so in equation 6n the ratio is limited to be .

The heat flux into the sheath is controlled by equation 6o [30]:

(6o) |

where here the sheath heat transmission coefficient is taken to be . Zero gradient (Neumann) boundary condition is used for the electron temperature, so the loss of energy at the sheath is implemented as an additional removal of thermal energy from the final grid cell.

The boundary condition on the electron density at the sheath entrance should be as unrestrictive as possible, to avoid over-constraining the system of equations. Free boundary conditions have been tried, in which the second derivative at the boundary is zero, but zero gradient boundary conditions have been found to be more robust and so are used here.

### 2.3 Coordinate system

A Clebsch coordinate system is used, aligned with the magnetic field such that the equilibrium magnetic field is given by . The coordinate is a flux coordinate in the radial direction, is an angular coordinate, and the coordinate is aligned with the equilibrium magnetic field. In BOUT++ the metric tensor is assumed to be constant in the direction. Two orientations of this coordinate system are used here: For simulations of the large aspect-ratio ISTTOK device in section 4 the direction is the poloidal angle. For tokamak simulations in section 5 the standard BOUT/BOUT++ coordinates are used, in which is the toroidal angle. For long time simulations the large scale magnetic field might be expected to evolve, but this is not currently allowed for in this model: Parallel derivatives are taken along the starting magnetic field, and the deviation due to the electromagnetic potential is assumed to be small.

## 3 Numerical methods

The numerical schemes used to solve the model equations have been chosen for their conservation properties. Some care must be taken over the form of the equations used, since not all analytically equivalent forms have the same numerical properties. Particular care must be taken over the form and choice of numerical method for energy exchange terms, which convert one form of energy to another [17]. The shear Alfvén wave dynamics along the magnetic field includes exchange of energy between electromagnetic and energy, so that in the terms

(6p) |

the operators and should numerically obey the identity . Fortunately this is quite easy to achieve, and the standard central differencing schemes have this property. A similar relation exists for the sound wave coupling between parallel ion flow and electron pressure, with the same solution.

The exchange of energy between energy and electron thermal energy appears as divergence of diamagnetic current in the vorticity equation, and compression of flow in the pressure equation:

(6qa) | |||||

(6qb) | |||||

(6qc) |

Here the second form of the compression term is used (equation 6qc) and both terms discretised with central differences since this then ensures that these terms combine into a divergence without relying on the properties of , for example its divergence in the numerical scheme.

The advection terms are discretised using the scheme illustrated in figure 1. The potential is first interpolated onto cell corners, in order to preserve a divergence-free flow in the absence of magnetic curvature. The velocity on each cell boundary is then calculated by taking the derivative of the potential along the boundary. Note that this (along with the interpolation) makes the method at best 2-order accurate.

The value of the field at either side of the boundary is determined by the choice of numerical scheme. This is done in each dimension independently, so each scheme must construct the values of at the left () and right () boundaries. The finite volume schemes implemented to calculate these values include 2-order Fromm and 4-order XPPM methods [36]. It has been found that use of the XPPM method results in a slow convergence of the implicit time integration scheme usually used (CVODE from the SUNDIALS suite [37]), and so here the Fromm method is used for advection by and magnetic drifts.

### 3.1 Poloidal flows

Poloidal flows due to drift are important for the Geodesic Acoustic Mode (GAM) oscillation [23, 38], and can play a role in the experimentally observed edge asymmetries and flow patterns [39, 40, 41]. In the tokamak coordinate system used here, the poloidal coordinate is transformed into the coordinate parallel to the magnetic field. The usual way to represent poloidal flows is to split the flow into two pieces:

(6qr) |

The first term is then represented as a Poisson bracket, whilst the compressional terms leading to GAM oscillations are contained in the second term. The advantages of this approach are that each of these terms can be discretised using the Arakawa method [42] which has minimal dissipation and respects the symmetries of the underlying equations. Unfortunately in general geometry it is extremely difficult to ensure that the resulting numerical method conserves particles since these terms must combine so as to obey the analytic integral relation. Here we adopt the approach also used in [10], and write the advection in flux-conservative form:

(6qsa) | |||||

(6qsb) | |||||

(6qsc) |

The first two terms describe the drift-plane motion, which is calculated using the scheme illustrated in figure 1. The third term describes radial flow due to poloidal electric fields, whilst the fourth (last) term describes poloidal flow due to radial electric fields. These poloidal flows are more difficult to implement than the drift-plane motion, particularly in the vicinity of an X-point where the cell corner is shared between eight cells rather than the usual four. Here we adopt a simpler scheme for the poloidal flows, in which the gradients of the potential are first calculated at cell centers, then interpolated to the cell boundaries to calculate the flux.

### 3.2 Numerical dissipation

The numerical methods employed here use collocated central differencing, and so some form of numerical dissipation is necessary to control zig-zag/chequerboard modes on the grid scale. This numerical dissipation must be carefully chosen so as not to introduce unphysical modes or instabilities. Perpendicular to the magnetic field we use classical and anomalous perpendicular diffusion of density, pressure, and vorticity (, and in equations 1-5). In the direction parallel to the magnetic field we use a combination of 4-order dissipation operators, corresponding to hyperviscosity terms in the flow variables (, ), and Added Dissipation [43] in the scalar variables (, , ). Added Dissipation is implemented as a flow between cells which is driven by the third derivative of the plasma pressure at the cell boundary. This conserves the flux of the advected quantity (density, thermal energy, momentum) whilst suppressing grid-scale oscillations due to the collocated numerical scheme.

## 4 Poloidal limited tokamak geometry

We first simulate turbulence in the ISTTOK device [24], a large aspect-ratio, poloidally limited tokamak. Major radius is cm, minor radius cm, toroidal magnetic field T, and safety factor . In order to accomodate the poloidal limiter, we use a coordinate system in which the direction is aligned with the poloidal direction , shifted so that the toroidal angle becomes identified with the direction along the magnetic field, :

(6qst) |

This coordinate system is illustrated in figure 2. In the core region the magnetic flux-surfaces are closed, so where the domain connects onto itself a twist-shift boundary using FFTs in the direction is used to map one end of the flux-tube to the other. This method, and the neglect of metric tensor variation in the poloidal () direction are possible only due to the large aspect ratio, so this coordinate system cannot yet be used in realistic X-point geometry (section 5.4).

In order to maintain a quasi-steady state profile, sources of heat and particles are needed close to the inner boundary. To achieve a specified core density and temperature, we use a Proportional-Integral (PI) feedback controller on the heating and particle sources in the core region. These sources are poloidally uniform, and are limited so that they can only be positive, preventing unphysical removal of heat or particles.

A simulation has been run with 68 radial points, 16 toroidal, and 256 poloidal points. This corresponds to a resolution of mm in the radial direction, and mm in the poloidal direction. The typical turbulence length scale is a multiple of mm at , so these simulations are probably not fully resolved. Higher resolution simulations will be required in order to carry out a quantitative validation exercise.

The plasma profiles develop along with the turbulence, and are not prescribed. A snapshot of the electron pressure is shown in figure a, showing a poloidal asymmetry. The pressure gradient gives rise to a Pfirsch-Schlüter current, which can be seen in figure b. Both pressure and parallel current profiles contain large fluctuations, of the same order as the time averaged value, underlining the importance of treating the background and fluctuations on the same footing in the plasma edge.

Fluctuations in plasma pressure as a function of time are shown in figure 4, in which the input power was increased at around ms and ms. Differences between the fluctuations in the closed field-line region inside the Last Closed Flux Surface (LCFS), and those outside the LCFS can be seen: Inside the LCFS (figure a) fluctuations are of similar amplitude on the inboard and outboard side, whilst outside the LCFS (figure b) fluctuations are larger on the outboard side, and a clear difference can be seen in the average pressure between the inboard and outboard side. This difference in average pressure between the inboard and outboard sides can be maintained in ISTTOK partly because these regions are not connected by parallel transport outside the LCFS, since field-lines intersect the poloidal limiter after travelling of a poloidal circuit.

More detailed study with a higher resolution grid, and comparison against experiment, will be the subject of future work.

## 5 Tokamak X-point geometry

Hermes allows simulation of both axisymmetric transport (section 5.2) and 3D turbulence (section 5.4) in tokamak X-point geometry. Here we describe the simulation procedure, and some features of the results. As with the ISTTOK asimulations, more detailed analysis is left to a more specific future publication.

The resolution of these simulations is () points in the (radial, parallel, toroidal) directions. Due to the field-aligned coordinate system, the effective poloidal resolution is much higher than this, being determined by the pitch of the magnetic field lines [18].

### 5.1 Coordinates

For X-point tokamak simulations the standard BOUT/BOUT++ coordinates are used. In terms of orthogonal toroidal coordinates these are [44, 18]

(6qsu) |

where is the poloidal arc length per radian (minor radius for circular cross-section), the major radius, the toroidal magnetic field and the poloidal magnetic field. As described in [18], the shifted metric method [45, 46] is used to reduce cell deformation due to magnetic shear.

### 5.2 Fluid transport

As discussed in section 1, 3D plasma turbulence on transport timescales is challenging. This can be made more difficult by the simulation starting conditions, which may be far from an equilibrium solution to the evolving equations 1-5, resulting in transient axisymmetric oscillations which are time consuming to evolve through. In order to reach a quasi-steady state more quickly, we first evolve only the 2D axisymmetric transport equations, without plasma currents or electric fields, using imposed cross-field anomalous transport coefficients in the spirit of 2D transport codes such as SOLPS [25], UEDGE [39, 47] or EDGE2D [40].

Using spatially uniform anomalous diffusion coefficients m/s and m/s the result is shown in figure 5.

Since there is no recycling of plasma at the target plates, the density falls from midplane to divertor, and is lower in the divertor than would be typical experimentally. The inclusion of neutral gas and plasma recycling is beyond the scope of this paper, but will be reported elsewhere. As with the ISTTOK simulations, a PI feedback controller is used to control particle and power source, so as to achieve the required core density and temperature.

### 5.3 Solution to in X-point geometry

Once the fluid transport simulation has reached quasi-steady state, we then allow evolution of the vorticity equation and parallel Ohm’s law, still in 2D (axisymmetric) mode.

The calculation of electric potential in this model requires inverting the elliptic equation 6g, which is a 2D problem on a curved surface embedded in the 3D domain. The operator is written in divergence form as

(6qsv) |

The metric components which couple the toroidal () and parallel projection of the poloidal component () are non-zero, so this is a 3D operator in the field-aligned coordinates used here (equation 6qsu).

The technique used in most previous BOUT [44] and BOUT++ [18] simulations was to neglect derivatives along the magnetic field, being small relative to the and derivatives in the ordering used (), and Fourier transform in the toroidal direction, this being possible when the Boussinesq approximation is used. This reduces the problem to a set of 1D equations in radial coordinate , which are then solved using the direct Thomas algorithm for tridiagonal systems. This method is computationally efficient, but for the mode the magnitude of the diagonal elements in this tridiagonal matrix becomes equal to the sum of the magnitudes of the off-diagonal elements, and the matrix is not strictly diagonally dominant. By comparing this solver with an iterative method, it was found that despite the relative error in the component on a single solve being of the order of , the iterative solver did not result in the growth of a numerical instability in cases where the direct solver did. By varying the tolerances in the iterative solver, it has been verified that this is not due to an effective smoothing error in the iterative solver. We therefore conclude that the direct solver is responsible, and should not be used for the component.

In addition to the issues discussed above, a further correction has been identified for simulations in X-point geometry: The flute ordering () used to justify dropping parallel derivatives is not sufficient close to the X-point, particularly for low toroidal mode numbers. Neglecting parallel (poloidal) derivatives does not on its own cause numerical instability, since it does not result in a spurious source of energy, but does significantly alter the solution, and introduces sharp gradients which can in turn lead to numerical problems in other terms. This is illustrated in figure 6, which shows the potential at a single time early in the development of the axisymmetric electrostatic potential in X-point geometry. Neglecting parallel derivatives in this coordinate system decouples grid points in the poloidal direction, in particular the coordinate lines passing either side of the X-point. Close to the X-point (marked with a box), this produces a discontinuity in the poloidal direction. This discontinuity is unphysical, and is resolved by retaining the derivatives in equation 6qsv.

This new solver does not yet handle finite toroidal mode numbers, so modes are inverted using the same solver as previously used in BOUT++, making the flute approximation to neglect parallel (poloidal) derivatives. Since this simplification is unlikely to be accurate for low modes in X-point geometry, here we retain only and modes, removing by simulating only of the torus. Note that this decomposition in toroidal modes is only possible here because we have made the Boussinesq approximation in equation 6g. If this approximation is not used then a 3D solution of becomes necessary in this coordinate system. An efficient means to do this is currently under development, and will be reported elsewhere.

Using this new axisymmetric field solver, the resulting evolution of the electrostatic potential in the core region is shown for two poloidal locations on the same flux surface in figure 7.

Three oscillation frequencies are apparent: (1) An oscillation during the first s with a frequency of approximately kHz, in which the potentials at the top and bottom of the plasma are out of phase; (2) A strongly damped oscillation with frequency around kHz during the first s, in which the potentials on the same flux surface are approximately in phase; and (3) a much slower oscillation with frequency around kHz, with potentials in phase. The simple analytical large aspect-ratio estimates of local wave frequencies are: shear Alfvén wave kHz; GAM frequency kHz; and the parallel sound wave kHz.

The electrostatic potential after ms on the time axis used in figure 7 (where currents are turned on at ) is shown in figure a. The electrostatic potential is relatively constant in the closed field-line region, as shown in the plot of radial electric field in figure 9. This quasi-steady state is associated with a parallel current shown in figure b. The Pfirsch-Schlüter current pattern is seen in the closed field-line region, balancing the current due to the magnetic drift (diamagnetic current divergence). In the SOL and Private Flux (PF) region parallel currents are seen to flow into the sheath in steady state, requiring closing flows through the vessel walls. This is because in the boundary conditions (equation 6n) we have assumed conducting walls.

Also shown in figure 9 is the radial electric field for a case with reversed toroidal field . This reverses the sign of the vertical magnetic drift , leading to a modification of the radial electric field in the plasma edge. In the SOL the electric field is only slightly modified, since the sheath strongly constrains the electric field.

### 5.4 Turbulence in X-point geometry

Finally, once the axisymmetric electric field has reached quasi-steady state we turn off the imposed anomalous transport (, ), extend the mesh in the toroidal domain to create a 3D simulation, and add small-scale random noise to the vorticity to seed instabilities. This is then run until it reaches a saturated turbulent state.

A timeseries of density fluctuations taken from the grid cell just outside the separatrix (which is located on a cell boundary) at the outboard midplane is shown in figure 10. The density shows an initial drop, as the imposed anomalous diffusion has been turned off, and the plasma fluctuations start at small amplitude and so provide little transport initially.

Averaging over the toroidal angle and time, the Root Mean Square (RMS) density fluctuations are calculated from the non-axisymmetric components of the density fluctuations (equation 6qsw).

(6qsw) |

where represents an average over toroidal angle , and an average over toroidal angle and time. This is shown in figure 11, showing a ballooning character to the fluctuations, but with significant fluctuations at all poloidal angles close to the separatrix.

As can be more clearly seen in figure b, relatively large fluctuations are observed in the private flux region, in particular in the bad curvature region of the inner leg, close to the divertor plate. This feature has been predicted in gradient-driven simulations [48] and observed experimentally [49]. The drive for these fluctuations is probably a combination of interchange drive due to bad curvature in this region, and a conducting wall mode due to the sheath [50], though more detailed analysis is needed to confirm this.

Further analysis of longer timeseries will be the subject of future work. Here we have focussed on testing the ability of the Hermes model to simulate turbulent transport together with evolving profiles and electric fields. Since the anomalous transport coefficients chosen initially are unlikely to result in the same transport as the turbulence, the plasma profiles will continue to evolve. Long turbulent simulations are therefore still required to find a self-consistent solution, but the procedure described here greatly reduces the computational cost: The turbulence simulations (figure 10) required approximately core-hours (on Archer) per millisecond of simulated time, whereas the fluid simulations (figure 5) required approximately 16 core-hours per millisecond.

## 6 Conclusions

We have described the key features of Hermes, a new model based on BOUT++, which is being developed to study transport in the edge of magnetically confined plasmas, in particular tokamaks. Progress has been made towards self-consistent simulation of turbulence and profiles on transport timescales: a drift-reduced model has been chosen and analysed for its conservation properties; numerical methods have been developed which conserve integral properties of the analytical model and so allow long-time simulation of plasma oscillations; and a new method for solving the electrostatic potential in X-point geometry with a field-aligned coordinate system has been developed and implemented. Together, these improvements enable simulations to be carried out in poloidal limiter (ISTTOK) and diverted tokamak configurations, which self-consistently evolve the large-scale electric fields and currents alongside the turbulence. Significant work remains to be done, some of which has been discussed in previous sections. Verification and validation of such a complex model will take significant effort, and is ongoing. As part of this work, we plan to carry out comparisons against ISTTOK in the near term as part of a EUROfusion project. The description of transport in the tokamak edge is strongly influenced by interaction with neutral gas, so in parallel with the work described here we have developed neutral gas models which will be described elsewhere. Hot ion effects, primarily ion diamagnetic drift and parallel viscosity, will modify the results presented here. Including these effects introduces complications, particularly in the vorticity equation, but is currently under developement. Finally, the description of the radial electric field (poloidal flow) in tokamak plasmas is a complex and subtle topic, and the present treatment will need further refinement.

## Acknowledgements

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The authors gratefully acknowledge the support of the UK Engineering and Physical Sciences Research Council (EPSRC) under grant EP/K006940/1, and Archer computing resources under Plasma HEC consortium grant EP/L000237/1.

## References

### Footnotes

- : Plasma Phys. Control. Fusion

### References

- A S Kukushkin, H D Pacher, V Kotov, G W Pacher, and D Reiter. Fusion Engineering and Design, 86:2865–2873, 2011.
- F Romanelli et al. Fusion Electricity - A roadmap to the realisation of fusion energy. Technical report, EFDA, 2012.
- R Wenninger, F Arbeiter, J Aubert, L Aho-Mantila, R Albanese, R Ambrosino, C Angioni, J.-F Artaud, M Bernert, E Fable, A Fasoli, G Federici, J Garcia, G Giruzzi, F Jenko, P Maget, M Mattei, F Maviglia, E Poli, G Ramogida, C Reux, M Schneider, B Sieglin, F Villone, M Wischmeier, and H Zohm. 55:063003, 2015.
- V Naulin. J. Nucl. Materials, 363-365:24–31, 2007.
- Ph Ghendrih, C Norscini, F Hasenbeck, G Dif-Pradalier, J Abiteboul, T Cartier-Michaud, X Garbet, V Grandgirard, Y Marandet, and Y Sarazin. J.Phys.: Conf. Ser., 401:012007, 2012.
- D A D’Ippolito, J R Myra, and S J Zweben. Physics of Plasmas, 18:060501, 2011.
- P Ricci, F D Halpern, S Jolliet, J Loizu, A Mosetto, A Fasoli, I Furno, and C Theiler. Plasma Phys. Control. Fusion, 54:124047, 2012.
- F D Halpern, P Ricci, S Jolliet, J Loizu, J Morales, A Mosetto, F Musil, F Riva, T M Tran, and C Wersal. J. Comput. Phys., 315:388–408, 2016.
- P Tamain, Ph Ghendrih, E Tsitrone, V Grandgirard, X Garbet, Y Sarazin, E Serre, G Ciraolo, and G Chiavassa. J. Comput. Phys., 229(2):361–378, 2010.
- P Tamain, H Bufferand, G Ciraolo, C Colin, D Galassi, Ph. Ghendrih, F Schwander, and E Serre. J. Comput. Phys., In Press, 2016.
- P J Catto, A N Simakov, F I Parra, and G Kagan. Plasma Phys. Control. Fusion, 50:115006, 2008.
- A B Mikhailovskii and V S Tsypin. Beitr. Plasmaphys., 24:335â354, 1984.
- D Pfirsch and D Correa-Restrepo. Plasma Phys. Control. Fusion, 38:71–101, 1996.
- A B Mikhailovskii et al. Plasma Phys. Rep., 23:844, 1997.
- D Reiser. Physics of Plasmas, 19:072317, 2012.
- A N Simakov and P J Catto. Physics of Plasmas, 10(12):pp. 4744–4757, December 2003.
- B Scott. Physics of Plasmas, 10:963, 2003.
- B D Dudson et al. Comp. Phys. Comm., 180:1467–1480, 2009.
- B D Dudson et al. 81(01):365810104, 2015. doi:10.1017/S0022377814000816.
- Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano Zampini, Hong Zhang, and Hong Zhang. PETSc Web page. http://www.mcs.anl.gov/petsc, 2016.
- Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschelman, Lisandro Dalcin, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, Lois Curfman McInnes, Karl Rupp, Barry F. Smith, Stefano Zampini, Hong Zhang, and Hong Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.7, Argonne National Laboratory, 2016.
- Satish Balay, William D. Gropp, Lois Curfman McInnes, and Barry F. Smith. Efficient management of parallelism in object oriented numerical software libraries. In E. Arge, A. M. Bruaset, and H. P. Langtangen, editors, Modern Software Tools in Scientific Computing, pages 163–202. Birkhäuser Press, 1997.
- N Winsor, J L Johnson, and J M Dawson. Phys. Fluids, 11:2448, 1968.
- C Silva, I Nedzelskiy, H Figueireda, R M O Galvao, J A C Cabral, and C A F Varandas. Nucl. Fusion, 44:799–810, 2004.
- R Schneider, X Bonnin, K Borrass, D P Coster, H Kastelewicz, D Rieter, V A Rozhansky, and B J Braams. Contrib. Plasma Phys., 46(1-2):3–191, 2006.
- J R Angus and M V Umansky. Physics of Plasmas, 21:012514, 2014.
- J T Omotani, F Militello, L Easy, and N Walkden. Plasma Phys. Control. Fusion, 58:014030, 2016.
- B Scott. Physics of Plasmas, 12:102307, 2005.
- R Chodura. Phys. Fluids, 25:1628, 1982.
- P C Stangeby. The Plasma Boundary of Magnetic Fusion Devices. IoP, 2000.
- J Loizu, P Ricci, F D Halpern, and S Jolliet. Physics of Plasmas, 19:122307, 2012.
- M U Siddiqui, D S Thompson, C D Jackson, J F Kim, N Hershkowitz, and E E Scime. Physics of Plasmas, 23:057101, 2016.
- S Togo, T Takizuka, M Nakamura, K Hoshino, K Ibano, T Long Lang, and Y Ogawa. J. Comput. Phys., 310:109–126, 2016.
- Ph Ghendrih, K Bodi, H Bufferand, G Chiavassa, G Ciraolo, N Fedorczak, L Isoardi, A Paredes, Y Sarazin, and E Serre. Plasma Phys. Control. Fusion, 53(5):054019, 2011.
- H Bufferand, G Ciraolo, G Dif-Pradalier, P Ghendrih, Ph Tamain, Y Marandet, and E Serre. Plasma Phys. Control. Fusion, 56:122001, 2014.
- J L Peterson and G W Hammett. SIAM J. Sci. Computing, 35(B576), 2013.
- A C Hindmarsh et al. ACM Transactions on Mathematical Software, 31(3):363–396, 2005.
- D Zhou. Physics of Plasmas, 22:092504, 2015.
- T D Rognlien, D D Ryutov, N Mattor, and G D Porter. Physics of Plasmas, 6:1851, 1999.
- A V Chankin, J P Coad, G Corrigan, S J Davies, S K Erents, H Y Guo, G F Matthews, G J Radford, J Spence, P C Stangeby, and A Taroni. Contrib. Plasma Phys., 40(3-4):288–294, 2000.
- A V Chankin, G Corrigan, M Groth, P C Stangeby, and JET contributors. Plasma Phys. Control. Fusion, 57:095002, 2015.
- A Arakawa. J. Comput. Phys., 1(1):119–143, 1966.
- J Y Murthy. Numerical methods in heat, mass and momentum transfer. Lecture notes ME608, 2002, Purdue University.
- X Q Xu, M V Umansky, B Dudson, and P B Snyder. Comm. in Comput. Phys., 4(5):pp. 949–979, November 2008.
- A M Dimits. Phys. Rev. E, 48(5):4070–4079, Nov 1993.
- B Scott. Physics of Plasmas, 8(2):447, 2001.
- T D Rognlien, X Q Xu, and A C Hindmarsh. J. Comput. Phys., 175:249–268, 2002.
- M V Umansky, T D Rognlien, and X Q Xu. J. Nucl. Materials, 337-339:266–270, 2005.
- J R Harrison, G M Fishpool, and B D Dudson. J. Nucl. Materials, 463:757–760, 2015.
- D D Ryutov and R H Cohen. Contrib. Plasma Phys., 44(1-3):168â175, 2004.