2 The codes EUTERPE and ORB5

Both EUTERPE (Jost et al. Reference Jost, Tran, Appert, Cooper and Villard2000) and ORB5 (Tran et al. Reference Tran, Appert, Fivaz, Jost, Vaclavik and Villard1999) share a common origin. They were initially developed at SPC Lausanne. Both are PIC $\mathrm {delta}-f$ codes, and many features first developed for one of them were afterwards implemented into the other one, as a result of the close collaboration of the corresponding development groups at Greifswald and SPC/Garching, respectively. However, they also have significant differences, the most important one of which is that ORB5 handles only axisymmetric equilibrium, whereas EUTERPE was designed from the beginning for stellarators and can treat any magnetohydrodynamic (MHD) equilibrium calculated with VMEC, either tokamak or stellarator. Both codes also differ in the equilibrium distribution function that they use. The EUTERPE uses a local Maxwellian, while ORB5 can use either a local or a canonical Maxwellian (Angelino et al. Reference Angelino, Bottino, Hatzky, Jolliet, Sauter, Tran and Villard2006; Vernay et al. Reference Vernay, Brunner, Villard, McMillan, Jolliet, Tran, Bottino and Graves2010). The ORB5 is also more flexible in the scheme used for the time evolution of the distribution function, it being able to use a direct $\mathrm {delta}-f$ scheme (Allfrey & Hatzky Reference Allfrey and Hatzky2003) or a standard $\mathrm {delta}-f$ one, while in EUTERPE only the standard $\mathrm {delta}-f$ scheme has been implemented so far. More details on the current status of the ORB5 code can be found in Lanti et al. (Reference Lanti, Ohana, Tronko, Hayward-Schneider, Bottino, McMillan, Mishchenko, Scheinberg, Biancalani, Angelino and Brunner2019).

As a matter of fact, the applications of ORB5, restricted to the tokamak domain, have been more extensive than those of EUTERPE, and many nonlinear turbulence tokamak simulations with ORB5 have been reported (see for instance Villard et al. (Reference Villard, McMillan, Lanti, Ohana, Bottino, Biancalani, Novikau, Brunner, Sauter and Tronko2019) and references therein). In the case of EUTERPE, nonlinear electrostatic simulations have only been studied in a screw pinch geometry so far (Sánchez et al. Reference Sánchez, Kleiber, Hatzky, Soba, Sáez, Castejon and Cela2010).

Both codes solve the gyrokinetic equations derived from a variational principle (Hahm Reference Hahm1988; Tronko, Bottino & Sonnendrucker Reference Tronko, Bottino and Sonnendrucker2016). In this work, we restrict ourselves to the electrostatic case, although both codes have implemented equations for the evolution of the electromagnetic potential.

An equation for the evolution of the distribution function of each kinetic species, $a$, is solved:

(2.1)\begin{equation} \frac{\partial f_a}{\partial t}+ \dot{\boldsymbol{R}} \frac{\partial f_a}{\partial \boldsymbol{R}} + \dot{v}_{\|} \frac{\partial f_a}{\partial v_{\|}} = \sum_b {\mathsf {C}}(\,f_a, f_b) + S, \end{equation}

where $S$ here represents the sources that can be added to control the numerical noise and/or sustaining density and temperature profiles; $S=S_{\textrm {ka}} + S_{\textrm {ha}}$, with $S_{\textrm {ka}}$ being a Krook operator and $S_{\textrm {ha}}$ a heating source term. A Krook operator, including corrections for energy, momentum and ZF conservation, is implemented in ORB5 (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008) while in EUTERPE only a simple Krook operator without corrections is implemented.

In this work, collisionless simulations (with $\sum _b {\mathsf {C}}(\,f_a, f_b)=0$) are analysed, except those in § 5, for which a pitch angle scattering collision operator is used (Kauffmann et al. Reference Kauffmann, Kleiber, Hatzky and Borchardt2010; García-Regaña et al. Reference García-Regaña, Kleiber, Beidler, Turkin, Maaßberg and Helander2013).

The equations of motion, $\dot {\boldsymbol {R}}$ and $\dot {v}_{\|}$, are

(2.2)\begin{equation} \left.\begin{array}{c} \dot{\boldsymbol{R}}=\underbrace{v_{\|} \boldsymbol{b} + \dfrac{\mu B + v_{\|}^{2}} {B^*_{a} \varOmega_a} \boldsymbol{b} \times \boldsymbol{\nabla} \boldsymbol{B} + \dfrac{v_{\|}^2 } {B^*_{a} \varOmega_a} (\boldsymbol{\nabla}\times \boldsymbol{B})_{\perp} - { \dfrac{ \boldsymbol{\nabla} \phi _{\textrm{ext}} \times \boldsymbol{b} } {B^*_{a}}} }_{{\dot{\boldsymbol{R}}^{(0)}}} {\underbrace{- \dfrac{ \boldsymbol{\nabla} \langle \phi \rangle \times \boldsymbol{b} }{B^*_{a}}}_{{\dot{\boldsymbol{R}}^{(1)}}}},\\ \dot{v}_{\|} = \underbrace{- {\mu} \left[ \boldsymbol{b} + \dfrac{v_{\|}} {B^*_{a} \varOmega_a}(\boldsymbol{\nabla} \times \boldsymbol{B})_{\perp} \right] \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B}}_{{\dot{v}_{\|}^{(0)}}} { \underbrace{- \dfrac{q_{a}} {m_{a}} \left[ \boldsymbol{b} + \dfrac{ v_{\|} } {B^*_{a} \varOmega_{a}} \left(\boldsymbol{b} \times \boldsymbol{\nabla} \boldsymbol{B} + (\boldsymbol{\nabla} \times \boldsymbol{B})_{\perp} \right)\right] \boldsymbol{\cdot} \boldsymbol{\nabla} \langle \phi \rangle }_{{\dot{v}_{\|}^{(1)}}}} , \end{array}\right\} \end{equation}

where $\boldsymbol {R}$ is the position of the gyrocentre and $v_{||}$ the parallel velocity; and the upper dot means time derivative. The magnetic moment, $\mu$, is a constant of motion ($\dot {\mu }=0$); $q_a$ and $m_a$ are the charge and mass, respectively, of the species $a$; $\varOmega _a={q_aB}/{m_a}$ and $B^*_{a}= B + ({m_a v_{||}}/{q_a}) \boldsymbol {b} \boldsymbol {\cdot } \boldsymbol {\nabla } \times \boldsymbol {b}$. The fourth term in the right-hand side of the equation for $\dot {\boldsymbol {R}}$ represents the contribution of a long-wavelength ambient electric field $\boldsymbol {\nabla }\phi _{\textrm {ext}}$ that can be introduced externally and is kept fixed during the simulation, while the fifth one represents that of the perturbed electrostatic potential $\phi$, which is consistently evolved.

A $\mathrm {delta}-f$ splitting is used in both codes so that the distribution function is separated into two parts,

(2.3)\begin{equation} f_a(\boldsymbol{R}, v_{||}, \mu, t) = f_{\textrm{Ma}}(\boldsymbol{R}, v_{||}, \mu)+ \delta f_a(\boldsymbol{R}, v_{||}, \mu, t), \end{equation}

with $f_{\textrm {Ma}}$ being a local Maxwellian in EUTERPE, and either a local or a canonical Maxwellian distribution function in ORB5.

With this splitting, a general nonlinear equation for $\delta f$ can be obtained,

(2.4)\begin{equation} \frac{\partial \delta f_a}{\partial t} + \dot{\boldsymbol{R}} \frac{\partial \delta f_a}{\partial \boldsymbol{R}} + \dot{v}_{\|} \frac{\partial \delta f_a}{\partial v_{\|}} = - { {\dot{\boldsymbol{R}}^{(1)}}} \frac{\partial f_{\textrm{Ma}}}{\partial \boldsymbol{R}} - {{\dot{v}_{\|}^{(1)}}} \frac{\partial f_{\textrm{Ma}}}{\partial v_{\|}} + \sum_b {\mathsf {C}}(\,f_a, f_b) + S, \end{equation}

assuming that $f_{\textrm {Ma}}$ is the equilibrium distribution function, or control variate (Aydemir Reference Aydemir1994). This equation can be linearized if $\dot { \boldsymbol {R}}^{(1)}$ and $\dot {v}_{\|}^{(1)}$ are dropped from $\dot { \boldsymbol {R}}$ and $\dot {v}_{\|}$ in the left-hand side of (2.4).

The distribution function is discretized by using a Klimontovich representation with markers, or quasiparticles, which follow the particle trajectories and carry a contribution to the distribution function (weight) whose evolution equation can be obtained from the equation for $\delta f$,

(2.5)\begin{equation} \frac{\textrm{d}w_p}{\textrm{d}t}=\left. \frac{\varOmega_p}{N_s} \frac{\textrm{d} \delta f_a}{\textrm{d} t} \right\vert_{Z_p} = \frac{\varOmega_p}{N_s} \left[- {{ \dot{ \boldsymbol{R}}^{(1)}}} \frac{\partial f_M}{\partial \boldsymbol{R}} - {{ \dot{v}_{\|}^{(1)}}} \frac{\partial f_M}{\partial v_{\|}} + \sum_b \textsf{\textit{C}}(\,f_a, f_b) + S \right]_{Z_p}, \end{equation}

where $\varOmega _p$ is the phase space volume associated with the p-th marker, $w_p$ is its weight, $N_s$ is the total number of markers in the simulation and $Z_p$ represents the location of the marker in phase space.

The system of equations is closed with the quasi-neutrality equation, which for the electrostatic case with adiabatic electrons, and using a long-wavelength approximation, reads

(2.6)\begin{equation} q_i \langle n_i \rangle - \frac{e n_0 (\phi - \{\phi\}_s)}{T_e}= - \boldsymbol{\nabla} \frac{m_in_0}{B^2} \boldsymbol{\nabla}_{\perp} \phi. \end{equation}

Here, the $\langle \rangle$ represents a gyro-average, and $\{ \}_s$ represents a flux-surface average.

In EUTERPE, cylindrical coordinates are used for pushing the particles (markers) and PEST (also known as $\theta ^*$) magnetic coordinates (Grimm, Greene & Johnson Reference Grimm, Greene and Johnson1976) are used for the fields. The ORB5 uses PEST magnetic coordinates both for the fields and also for pushing the particles.

For more details about the equations solved and their implementation in each code, the reader is referred to Jost et al. (Reference Jost, Tran, Cooper, Villard and Appert2001), Kornilov et al. (Reference Kornilov, Kleiber, Hatzky, Villard and Jost2004), Jolliet et al. (Reference Jolliet, Bottino, Angelino, Hatzky, Tran, Mcmillan, Sauter, Appert, Idomura and Villard2007), Slaby et al. (Reference Slaby, Könies, Kleiber and García-Regaña2018) and Lanti et al. (Reference Lanti, Ohana, Tronko, Hayward-Schneider, Bottino, McMillan, Mishchenko, Scheinberg, Biancalani, Angelino and Brunner2019).

3 Benchmark EUTERPE–ORB5

The code EUTERPE has been benchmarked in linear settings against several codes, such as TORB (Sánchez et al. Reference Sánchez, Kleiber, Hatzky, Soba, Sáez, Castejon and Cela2010), GENE (Helander et al. Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015; Görler et al. Reference Görler, Tronko, Hornsby, Bottino, Kleiber, Norscini, Grandgirard, Jenko and Sonnendrücker2016), GYSELA, GKW and ORB5 (Görler et al. Reference Görler, Tronko, Hornsby, Bottino, Kleiber, Norscini, Grandgirard, Jenko and Sonnendrücker2016) and XGC (Cole et al. Reference Cole, Hager, Moritaka, Lazerson, Kleiber, Ku and Chang2019). In a turbulence nonlinear setting it has only been compared with TORB for electrostatic simulations in a screw pinch geometry, so far (Sánchez et al. Reference Sánchez, Kleiber, Hatzky, Soba, Sáez, Castejon and Cela2010).

As a step previous to the nonlinear simulations in stellarator geometry with EUTERPE, we carry out a benchmark against ORB5 in tokamak geometry. We choose a tokamak equilibrium, matching the well known cyclone base case (CBC) tokamak equilibrium (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) at half-radius. The main parameters of this equilibrium are major radius $R=1.7\ \textrm {m}$, minor radius $a=0.625\ \textrm {m}$, the magnetic field at the axis $B(0)=1.91\ \textrm {T}$ and the safety factor $q(0)=0.85$ and $q(a)=3.24$.

A set of linear and nonlinear simulations were carried out with both codes using this MHD equilibrium, adiabatic electrons, and ideal density and temperature profiles constructed according to the analytic formula

(3.1)\begin{equation} X=X_* \exp \left\{\frac{\kappa_X \varDelta_X}{2} \ln \left[ \frac{\textrm{cosh}\left(\dfrac{r/a-(\rho_0+\varDelta_X)}{\varDelta_s}\right)}{\textrm{cosh}\left(\dfrac{r/a-(\rho_0-\varDelta_X)}{\varDelta_s}\right)}\right] \right\}, \end{equation}

where $X=\{n_i, n_e,T_e,T_i\}$ represents a radial profile of density or temperature for ions or electrons. Here $X_*$ is the value of $X$ at the reference normalized radius $\rho _0=0.5$, and $\kappa _X$ and $\varDelta _X$ are specific parameters for each profile. An example of this type of profile is shown in figure 1.

Figure 1.

Density and temperature profiles (a) and their characteristic scale lengths (b) for the analytical profiles defined in (3.1) with parameters $\kappa _{\textrm {Ti}} = \kappa _{\textrm {Te}}=3.774$, $\kappa _{\textrm {ni}}=\kappa _{\textrm {ne}}=0.807$, $\varDelta _n=\varDelta _{\textrm {Te}}=\varDelta _{\textrm {Ti}}=0.4$, $\varDelta _s=0.04$.

3.1 Linear simulations

With these model density and temperature profiles we carry out a set of linear simulations for several values of the parameter $\kappa _{\textrm {Ti}}=2,2.54,3,3.5$. In all cases, the electron temperature and density gradients are kept the same with $\kappa _{\textrm {Te}}=0.3775$, $\varDelta _s=0.04$, $\kappa _{\textrm {ni}}=\kappa _{\textrm {ne}}=0.8$, $\varDelta _{\textrm {ne}}=\varDelta _{\textrm {Te}}=0.4$ and $\varDelta _{\textrm {ni}}=\varDelta _{\textrm {Ti}}=0.2$. Simulations are global in radius, and the maximum instability appears around $r/a\sim 0.38$. Note that the growth rate of ion temperature gradient (ITG) modes depends on the local value of the ITG, which is maximum around $r/a=0.12$, but the modes have a radial extent and consequently are affected by the profiles and magnetic structure in the radial region they cover. The same settings are used in ORB5 and EUTERPE codes: a spatial resolution in the radial,Footnote ¹$s$, poloidal, $\theta ^*$, and toroidal, $\phi$, directions $n_s \times n_{\theta ^*} \times n_{\phi }=64 \times 128\times 64$ and a squared low pass Fourier filter keeping modes $-63 < m <63$, plus a diagonal filter suppressing modes with $|m - n/{\raise.3pt-\kern-5pt\iota}|> 5$. Only a toroidal mode $n = 16$ is kept in both codes to simplify the comparison. This mode is the most unstable one when a simulation including several $n$ modes is carried out.

The growth rate and the real frequency are obtained by fitting the time evolution of the potential to an oscillation whose amplitude grows in time. From the fitting, both the growth rate and real frequency of the most unstable mode can be obtained. The results of these fittings for both codes are shown in figure 2. As shown in the figure, the agreement with the growth rates is excellent, while it is not that good in the real frequencies, which are more sensitive to the precise time evolution of the modes. Note that in these simulations several $m$ modes are present with comparable growth rates and different frequencies. This can affect the fitting to a unique frequency. Also, differences in the code implementation details can affect the frequency results, while the growth rates appear to be more robust.

Figure 2.

Comparison of growth rate (a) and frequency (b) of the most unstable mode in ORB5 and EUTERPE for different values of $\kappa _{\textrm {Ti}}$.

3.2 Nonlinear simulations

Once the results of both codes have been successfully compared for linear simulations we proceed with the comparison in a nonlinear setting. We use the same CBC magnetic equilibrium, and for the density and temperature profiles we use analytic profiles described by (3.1) with parameters from (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008). The density and temperature profiles are the same for electrons and ions with $\kappa _{\textrm {Ti}} = \kappa _{\textrm {Te}}=3.774$, $\kappa _{\textrm {ne}}=\kappa _{\textrm {ni}}=0.807$ and $\varDelta _n=\varDelta _{\textrm {Te}}=\varDelta _{\textrm {Ti}}=0.04$. These profiles are shown together with their characteristic scale lengths in figure 1.

For this comparison we target the full-volume-integrated heat flux, $Q_i$, as the fundamental physical quantity to study. This quantity is obtained by averaging the radial flux of kinetic energy over all the markers used in the simulation,

(3.2)\begin{equation} Q_{i}=\sum_{j=1}^{N}(\,f_{0j}+w_j)\frac{m_j v_j^2}{2} \frac{\left\langle{\boldsymbol{E}} \right\rangle\times {\boldsymbol{B}}}{BB^*_{\|}}\cdot \frac{\boldsymbol{\nabla} s} {|\boldsymbol{\nabla} s|}. \end{equation}

Here $f_{0j}$ represents the average of the equilibrium distribution function on the $j\text {-}\mathrm {th}$ marker's volume and $w_j$ is the contribution to the $\delta f$ (weight) of the $j\text {-}\mathrm {th}$ marker.

In addition to the heat flux, the volume-averaged heat conductivity, which is computed as

(3.3)\begin{equation} \chi_i = \frac{\sum_k V_k \chi_{ik}}{ \sum_k V_k}, \end{equation}

is also compared between both codes. The computation of this magnitude requires a radial binning to obtain the conductivity in a finite number of radial positions, with $V_k$ being the volume of the $k-\mathrm {th}$ radial bin.

The values of the heat conductivity, $\chi _{ik}$, are obtained by averaging properties of markers within each radial bin, as $\chi _{\textrm {ik}} = Q_{\textrm {ik}} /\boldsymbol {\nabla } T_{\textrm {ik}}$, with $\boldsymbol {\nabla } T_{ik}$ being the average ITG in the $k\text {-}\mathrm {th}$ bin, and

(3.4)$$ \begin{equation} Q_{\textrm{ik}} = \sum_{s_{k_0} \leq s_j < s_{k_1} }(\,f_{0j}+w_j) \frac{m_j v_j^2}{2} \frac{\left\langle{\boldsymbol{E}}\right\rangle\times {\boldsymbol{B}}}{BB^*_{\|}}\cdot \frac{\boldsymbol{\nabla} s} {|\boldsymbol{\nabla} s|}, \end{equation}$$

where $s_j$ represents the radial location of the $j\text {-}\mathrm {th}$ marker, and $s_{k_0}$ and $s_{k_1}$ represent the limits of the $k\text {-}\mathrm {th}$ radial bin, centred around the radial position $s_k$.

In order to get a stable signal for the heat flux, we resort to the Krook operator which is now implemented in both codes and has proved to allow the stabilization of noise and heat flux in nonlinear simulations with ORB5 (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008). In EUTERPE only the simple Krook operator is implemented, while in ORB5 a corrected Krook, including energy and momentum conservation corrections, is also available. Then, for a fair comparison, we switch off these corrections in ORB5 to use the simple Krook version in it also. No additional sources or weight smoothing is used in these simulations. The same input parameters are used in both codes: a spatial resolution $n_s \times n_{\theta ^*} \times n_{\phi }=256 \times 256\times 128$, a squared low pass Fourier filter allowing modes $-63 < m <63$, $-63 < n <63$ and a diagonal filter suppressing modes with $|m - n / {\raise.3pt-\kern-5pt\iota} | > 5$; quadratic splines are used in both cases. The comparison of the heat fluxes and conductivities is shown in figures 3 and 4 for two simulations carried out with a different number of markers, $150 \times 10^{6}$ and $450 \times 10^{6}$, respectively. The heat conductivity is normalized to gyro-Bohm units $\chi _{\textrm {GB}}=\rho _i c_s / a$, with $\rho _i$ the ion Larmor radius, $c_s=\sqrt {T_e/m_i}$ the sound speed and $a$ the minor radius.

Figure 3.

Full-volume integrated heat flux from (3.2) (a) and volume-averaged heat conductivity from (3.3) (b) for two simulations carried out with ORB5 and EUTERPE using the same settings (see the text) and using $150\times 10^6$ markers.

Figure 4.

Full-volume integrated heat flux from (3.2) (a) and volume-averaged heat conductivity from (3.3) (b) obtained in simulations carried out with ORB5 and EUTERPE using the same settings (see the text) and using $450 \times 10^6$ markers. Results of a simulation carried out with EUTERPE using the same initial noise in the markers’ weights as in ORB5 is also shown in dotted line (labelled EUT-noise $1\times 10^{-7}$).

In the first case, with $150 \times 10^6$ markers, which is shown in figure 3, the statistics is not enough, and the numerical noise is significant. This is manifested in decaying heat flux signals and translates also into a mismatch between ORB5 and EUTERPE results, particularly in the heat conductivity, which is affected by both errors in the estimation of heat flux in the radial bins and also by the error in the temperature gradient estimates.

Significantly better agreement between codes is obtained with an increased number of markers (figure 4) as compared with the previous one (figure 3). The agreement in the latter case is very good in both the heat flux and conductivity, in spite of the implementation details, which suggests that the statistics of $450 \times 10^6$ markers is large enough in this case.

A clear difference between ORB5 and EUTERPE simulations can be observed in the initial stages in figure 4: the exponential growth due to the linear instability initiates first in the EUTERPE case. The difference is due to the fact that a larger amplitude initial noise was used in the EUTERPE case as compared with the ORB5 one, and also to the slightly different implementations in both codes. The results of a simulation carried out with EUTERPE using $450 \times 10^6$ markers and a smaller amplitude initial noise (exactly equal to the value used in ORB5) are also shown in the figure (dotted line). It is clear that the exponential growth of heat flux and conductivity starts later than in the case with larger amplitude noise, although an exact match of the growth between both codes is not obtained due to the different implementations. However, independently of the initial noise used and the different numerical details in both codes, the saturated levels in heat flux and conductivity, which are the physically meaningful quantities, are the same in both codes. Note that despite the increase in the markers’ statistics from figure 3 to figure 4 a significant level of burstiness is still present in figure 4, which suggests that it is not related to the numerical noise but to the presence of phenomena with significant time and space correlation lengths, as expected for a turbulent sate (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008).

4 Characterization of NCSTs

In this section, we test the use of the NCSTs and study their influence on the quality and physical results of the simulations.

A simple Krook operator (Krommes Reference Krommes1999; McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008) has been implemented in EUTERPE without energy or momentum correction. This simple Krook operator is one of the source terms present in (2.4) with the form $S_{\textrm {ka}} = -\gamma _{\textrm {ka}} \delta f_a$, where $\gamma _{\textrm {ka}}$ represents the strength of the source for species $a$. As in this work we present simulations with only kinetic ions, in the following sections we will drop the $a$ index related to the kinetic species and $\gamma _k$ will refer to ions. This kind of source term is commonly used to stabilize collisionless simulations. In collisional cases, its use is more questionable as it can distort the effect of a more realistic collision operator. In any case, as a common sense rule $\gamma _{k}$ should be chosen small as compared with the growth rate of the unstable modes, and smaller than the collision frequency in collisional simulations, in order to not affect significantly the results. In this section, we will show how the size of $\gamma _{k}$ affects the results in different simulations.

In addition to the Krook operator, a weight smoothing scheme (Sonnendrücker et al. Reference Sonnendrücker, Wacher, Hatzky and Kleiber2015) has also been implemented in EUTERPE. For this purpose, the markers are classified in a quad tree (QT) according to their coordinates in phase space and the weights of neighbour markers are mixed after a prescribed number of time steps, so that their variance is reduced. The weights of a couple of neighbouring markers, $w_1$ and $w_2$, are redefined at each smoothing step as

(4.1)\begin{gather} w_1^{i+1}=(1-\textrm{e}^{-{d^2}/{2\sigma}})w_1^i+\textrm{e}^{-{d^2}/{2\sigma}}\frac{w_1^i+w_2^i}{2} \end{gather}

(4.2)\begin{gather}w_2^{i+1}=(1-\textrm{e}^{-{d^2}/{2\sigma}})w_2^i+\textrm{e}^{-{d^2}/{2\sigma}}\frac{w_1^i+w_2^i}{2} , \end{gather}

where the $i$ super-index indicates quantities before the redefinition and $i+1$ the quantities after the mixing, and $d^2=(v_{||1}-v_{||2})^2 + (v_{\perp 1}-v_{\perp 2})^2$ represents the separation of markers in velocity space coordinates.

It is clear from the definition of smoothed weights in (4.1) and (4.2) that after the weights update, the density is conserved and the variance of the weights is reduced. However, with this smoothing scheme neither the energy nor the momentum of markers are conserved. Three parameters can be modified which affect the strength of the weight smoothing. First, the frequency at which the weight smoothing is applied, which we call $f_{\textrm {QT}}$. The second parameter is $\sigma$ in (4.1) and (4.2), which determines the strength of the smoothing and also affects the energy and momentum conservation. The smaller $\sigma$ is, the closer the markers have to be in order to mix their weights. As $\sigma$ increases, the exponential factor in (4.1) and (4.2) increases for large values of the markers separation $d$ and, consequently, the conservation of momentum and energy will be worse. In the following, we will call this parameter either $\sigma$ or $\sigma _{\textrm {QT}}$. Finally, there is a free parameter in the way the markers are sorted in the QT. It is the minimum number of markers that are required in each phase space cell. Note that cells have not the same width but are defined so that each of them contains at least a minimum number of markers per cell, which we call $N_{\textrm {QT}}$. The smoothing is only applied to neighbour markers within the same cell.

These tools have proved to help in mitigating the growth of numerical noise in ORB5, but they can also affect the ZF dynamics. While they have shown good results in ORB5 tokamak simulations they have not been previously tested in a global PIC code in stellarators.

In addition, a heating source term,

(4.3)\begin{equation} S_{\textrm{ha}} = -\gamma_{\textrm{ta}}(\delta f_a(s, \epsilon) -f_a^0(s, \epsilon)\frac{\tilde{n}_a(s)}{n_{a0}(s)})- \gamma_{\textrm{na}}f_a^0(s, \epsilon)\frac{\tilde{n}_a(s)}{n_{a0}(s)}, \end{equation}

is used, which allows the restoration of the density and temperature profiles on a long time scale, determined by the size of $\gamma _{\textrm {ta}}$ and $\gamma _{\textrm {na}}$ (McMillan et al. Reference McMillan, Jolliet, Tran, Villard, Bottino and Angelino2008). As in the case of the Krook operator, the time constants of this heating source term should be large enough as compared with the maximum growth rate and collision time, in order not to perturb the simulation results.

The quantities $f_a^0(s, \epsilon )$, $\delta f_a(s, \epsilon )$ represent, respectively, the average values of the equilibrium and perturbed distribution functions of species $a$ in a radial bin around the radial position $s$, and the kinetic energy $\epsilon$, of the marker. Here, $n_{a0(s)}$ and $\tilde {n}_a(s)$ are the equilibrium and perturbation density around this position. Note that $n_{a0}(s)=\sum _{\epsilon } f_a^0(s, \epsilon )$ and $\tilde {n}_a(s)=\sum _{\epsilon }\delta f_a(s, \epsilon )$. As for the Krook operator, we will drop the species index from $\gamma _{\textrm {ta}}$, $\gamma _{\textrm {na}}$ because we only deal with a kinetic species, the main ions.

4.1 Influence of NCSTs on linear growth rates in a tokamak

Before going to the characterization of NCSTs in stellarators, we first study how using these tools affects the growth rates of linearly unstable modes in tokamak geometry. We use the same CBC MHD equilibrium and analytical density and temperature profiles given in (3.1). First, we study the use of a Krook operator. We run a set of linear collisionless simulations using the same density and temperature profiles from figure 1 and using a Krook operator with different values of the constant $\gamma _k = 10^{-6}, 10^{-5}, 10^{-4}, 5\times 10^{-4}, 10^{-3}$, in units of $\varOmega _* ^{-1}$, with $\varOmega _* = e B_* / m$ the ion cyclotron frequency, and $B_*$ the average value of the magnetic field strength at the magnetic axis. The results of these simulations are shown in figure 5. The growth rate of the most unstable mode for the reference case, without Krook operator, is $\gamma _m=8.5 \times 10^{-4} \varOmega _*$. A shown in the figure, the Krook operator has a strong effect on the linear growth rate when $\gamma _k$ approaches the linear growth rate of the unstable modes. It can be inferred from these results that keeping $\gamma _k < \gamma _m /10$ does not alter the growth rate in more than 10 %.

Figure 5.

Linear growth rate of the most unstable mode in a set of linear simulations in the CBC MHD equilibrium using a Krook operator with different values of $\gamma _k$.

In figure 6 the influence of using the weight smoothing is shown. The parameter with a stronger influence on the growth rate is $\sigma _{\textrm {QT}}$, while the frequency and minimum number of markers per bin have a weak influence. This is easily understood because $\sigma _{\textrm {QT}}$ enters the expressions (4.1) and (4.2) in an exponential function determining the distance in phase space at which the smoothing is effective. Note the different scales used in the right-hand vertical axis.

Figure 6.

Influence of the weight smoothing on linear growth rate in the CBC equilibrium for different values of the frequency (a), minimum number of markers in each QT cell (b) and $\sigma _{\textrm {QT}}$ (c). Note the different scales used in the right-hand vertical axis. The reference parameters, used when only one of them is changed, are $N_{\textrm {QT}}=10$, $\sigma _{\textrm {QT}}=0.4$ and $f_{\textrm {QT}}=1$.

Finally, figure 7 shows the results of a set of linear simulations with the same parameters in which only the heating source term is used with different values of $\gamma _t$ and $\gamma _n$. We show the results obtained by changing each parameter separately. It is clear that the heating source term has a small influence on the linear growth rate even for values of $\gamma _t$ and $\gamma _n$ close to, or even slightly larger than, the maximum growth rate.

Figure 7.

Influence of the heating source term on the linear growth rate in the CBC equilibrium for different values of $\gamma _t$ (a) and $\gamma _n$ (b).

4.2 Influence of NCSTs on linear properties in a stellarator configuration

Now we turn to study the effect of the NCSTs in stellarator configurations. We will use the standard magnetic configuration of W7-X (Ref_11_EIM) and model density and temperature profiles given by

(4.4)\begin{equation} X=X_* \exp \left\{ \frac{-\kappa _x}{1-\cosh^2\left(\dfrac{\rho_0}{\varDelta _x}\right)} \left[\varDelta _x \tanh\left(\frac{\rho-\rho_0}{\varDelta _x}\right) - \cosh^2\left(\frac{\rho-\rho_0}{\varDelta _x}\right)\right] \right\}, \end{equation}

where $X=\{n, T_e, T_i\}$ represents the profiles for density, electron and ion temperatures, and $X_*$ their values at the reference position $\rho _0$, with $\rho =r/a$, $\kappa _{\textrm {Ti}} =3$, $\kappa _{n} =1$, $\rho _0=0.5$, $\varDelta _{\textrm {Ti}} = 0.1$, $\varDelta _n=0.3$. The electron temperature profile is flat with $T_e=T_*=11.5\ \textrm {keV}$ and $n_*=1\times 10^{19}\ \textrm {m}^{-3}$. The ion temperature profile has a maximum $T_i$ gradient at half-radius, $r/a=\rho _0=0.5$. For this equilibrium we run a linear collisionless simulation with adiabatic electrons, without any source term or weight smoothing, in which the growth rate of the most unstable mode is $\gamma _{\max }\sim 8.5\times 10^{-4}\varOmega _*$. This is used as a reference simulation without using any NCST.

4.2.1 Influence of NCSTs on the linear growth rate

In figure 8 the influence of Krook operator ($\gamma _k$), QT smoothing ($\sigma _{\textrm {QT}}$) and heating sources ($\gamma _t, \gamma _n$) on the linear growth rate of unstable modes is shown for a set of linear simulations in this W7-X configuration. We use the constant $\gamma _h$ for the cases with $\gamma _t=\gamma _n=\gamma _h$. As shown in the figure, and consistently with results from § 4.1, the effect of the Krook operator on the linear growth rate is important when $\gamma _k$ is close to the growth rate of the most unstable modes without any NCST used ($\gamma _{\max } = 8.5 \times 10^{-4} \varOmega _*$). The influence of the heating source is always smaller than 10 %, even for values of $\gamma _t$ and $\gamma _n$ similar or above the growth rate of the most unstable modes.

Figure 8.

Influence of the Krook operator (a), heating source, with equal values of $\gamma _t$ and $\gamma _n$, (b) and weight smoothing for $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$ (c) on the linear growth rate in a W7-X equilibrium (Ref_11_EIM) for different values of parameters.

With respect to the weight smoothing, figure 8 shows the increasing effect of weight smoothing on the linear growth rate with the parameter $\sigma _{\textrm {QT}}$. It is clear that for $\sigma _{\textrm {QT}}<0.2$ the effect on the growth rate is smaller than 10 %. There is a slight dependency of the effect on the growth rate with the minimum number of markers in the phase space cell, $N_{\textrm {QT}}$, with an increasing effect as this number is reduced (not shown). This is interpreted as a result of the sharp cutoff associated with the exponential term in (4.1) and (4.2). As the minimum number of markers in the bin is increased, the distance in phase space between markers whose weights are mixed is increased, and then the exponential factor decreases strongly and the effect of smoothing is then reduced. On the contrary, for reduced numbers of markers the distance between them is reduced and the smoothing becomes more important. As the smoothing should be done only between markers with similar velocity components in order to preserve, as much as possible, the energy and momentum conservation, a rule of thumb can be derived that the number $N_{\textrm {QT}}$ should be reduced to the minimum for better conservation properties.

No significant difference is found on the degree that the Krook operator and the heating sources affect the linear growth rate of unstable modes in tokamak and stellarator configurations. Concerning the weight smoothing, the reduction in the growth rate due to smoothing is larger in the stellarator configuration than in tokamak (see figures 6c and 8c). Part of the difference can be attributed to the different smoothing frequency, $\varOmega _*$, in both cases. Figure 8c corresponds to a larger smoothing frequency than that of figure 6, and then a larger influence of the smoothing can be expected. Besides, the reduction of growth rate with smoothing is larger in the stellarator case for the same smoothing frequency.

4.2.2 Influence of NCSTs on the linear evolution of ZFs

In this section, we study the influence of these tools on the linear evolution of ZFs in a stellarator configuration, both the residual level (Monreal et al. Reference Monreal, Calvo, Bustos and Kleiber2016) and the low-frequency characteristic oscillation (Mishchenko et al. Reference Mishchenko, Helander and Könies2008; Monreal et al. Reference Monreal, Sánchez, Calvo, Bustos, Parra, Mishchenko, Könies and Kleiber2017). We use the same standard configuration of W7-X as in § 4.2.1. In this case, we use flat density and temperature profiles with $T_i=T_e=5\ \textrm {keV}$. We run linear simulations of ZF relaxation for a very small radial wavelength of the ZF. These simulations are run with adiabatic electrons and are started with a perturbation to the density $\delta f$ that produces an initial potential $\phi \propto \cos (k_s \pi s)$, with $k_s\rho _i<0.1$. Under these conditions, the residual level is very small, in agreement with expectations (Monreal et al. Reference Monreal, Calvo, Bustos and Kleiber2016). First, we run a simulation with these parameters and without using any NCST, which we will use as a reference. In this simulation, we observe a low-frequency oscillation of the zonal potential with $\varOmega _{ZF} = 8 \times 10^{-6} \varOmega _*$. Then, we run a set of linear simulations with the same setting but using the NCSTs with different parameters. The results are shown in figure 9. The use of a Krook operator introduces a damping in the low-frequency ZF oscillations. Its effect increases with the value of $\gamma _k$, and it has an important effect even on the short-time evolution of the ZF oscillation for values $\gamma _k > \varOmega _{\textrm {ZF}}$. The effect of using a heating source on the ZF oscillation is smaller than that of the Krook operator, with a smaller effect on the damping of ZF oscillations, even for values of $\gamma _h$ slightly above the reference ZF frequency, $\gamma _h \sim \varOmega _{\textrm {ZF}}$. It should be noted that the ZF oscillation frequency depends on the magnetic geometry and is typically much smaller than the typical growth rate of unstable modes (Monreal et al. Reference Monreal, Sánchez, Calvo, Bustos, Parra, Mishchenko, Könies and Kleiber2017). In previous cases studied in W7-X, in § 4.2.1, the maximum linear growth rate was $\gamma _{\max } \sim 70\varOmega _{\textrm {ZF}}$ (at half-radius). Then, it can be expected that for a nonlinear simulation including realistic density, temperature profiles and ZF response, the parameter $\gamma _k$ of the Krook operator required for an effective noise control strongly affects the ZF evolution even at short times. The effect of heating sources can be expected to be slightly smaller, however.

Figure 9.

Influence of the Krook operator (a), heating source, with equal values of $\gamma _t$ and $\gamma _n$, (b) and weight smoothing (c) on the linear evolution of the zonal potential component in a W7-X equilibrium for different values of parameters $\gamma _k$, $\gamma _h$ ($\gamma _t=\gamma _n=\gamma _h$) and $\sigma _{\textrm {QT}}$. In (c) $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$.

The effect of weight smoothing on the ZF oscillations is very strong for all values of $\sigma _{\textrm {QT}}$ studied, from $0.1$ to $0.5$ (always with $N_{\textrm {QT}}=10, f_{\textrm {QT}}=1 \varOmega _*$). The long time linear evolution of ZFs depends on averages over the particle trajectories of quantities such as the radial drift or the radial displacement of trapped particles (see Monreal et al. Reference Monreal, Calvo, Bustos and Kleiber2016, Reference Monreal, Sánchez, Calvo, Bustos, Parra, Mishchenko, Könies and Kleiber2017). While the weight smoothing does not affect the trajectories of the markers, it changes the marker's weights, which has a cumulative effect that increases with time and can be expected to affect the long-time averages. However, it is interesting noting here that the effect on the short times is very small for any value of $\sigma _{\textrm {QT}}$ as compared with the effect of both Krook operator and the heating source term. The effect of using these tools on the short-time evolution of ZFs is shown in more detail in figure 10. It can be argued that the linear evolution of ZFs for long times as compared with the characteristic turbulence times (eddy turn-over time, or inverse of the growth rate) should not be very relevant in the turbulence saturation. In this perspective, the weight smoothing could be considered as a method of noise control with small effect on the ZF, while the influence on ZFs at short times of heating sources, and particularly that of the Krook operator, is larger, always depending on the $\gamma _k$ and $\gamma _h$ parameters values.

Figure 10.

Influence of the Krook operator (a), heating source, $\gamma _h$, (b) and weight smoothing, $\sigma _{\textrm {QT}}$ (with $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$) (c), on the linear evolution of the zonal potential component in W7-X equilibrium for different values of parameters.

A background long-wavelength radial electric field is known to have an effect on both the residual ZF level (Sugama & Watanabe Reference Sugama and Watanabe2009) and also on its oscillation frequency (Mishchenko & Kleiber Reference Mishchenko and Kleiber2012); then, it is worth checking the effect of using these NCST tools on the ZF evolution with a background electric field included in the simulation. The results of this test for two values of the background electric field are shown in figures 11 and 12. An electric field of the form $\phi '={\textrm {d}\phi }/{\textrm {d}s}=0.5$ is used in the first case, shown in figure 11, with $s=(r/a)^2$ being the radial coordinate, normalized toroidal flux. By comparing figures 8 and 11 we can appreciate the effect of the electric field on the ZF oscillations: it introduces a damping on the oscillation. For this value of the radial electric field, the change in this frequency and the residual ZF level is very small, however. The ZF oscillations are further damped when NCSTs are used, although the effect is slightly smaller than with no electric field (compare figure 11 and 8). With respect to the influence at short times, the same conclusion obtained without electric field can be extracted in this case: weight smoothing is much less distorting at short times, than using a Krook source term as well when an electric field is included.

Figure 11.

Influence of the Krook operator (a), heating source, $\gamma _h$, (b) and weight smoothing, $\sigma _{\textrm {QT}}$ (with $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$) (c), on the linear evolution of the zonal potential component in W7-X equilibrium for different values of parameters including an electric field ${\textrm {d}\phi }/{\textrm {d}r}=r/a^2$.

Figure 12.

Influence of the Krook operator (a), heating source, $\gamma _h$, (b) and weight smoothing, $\sigma _{\textrm {QT}}$ (with $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$) (c), on the linear evolution of the zonal potential component in a W7-X equilibrium for different values of parameters including an electric field ${\textrm {d}\phi }/{\textrm {d}r}=2r/a^2$.

In figure 12 we show the results for a series of simulations including a larger electric field, $\phi '={\textrm {d}\phi }/{\textrm {d}s}=1$. In this case, the electric field introduces a finite residual level, an increase in the oscillation frequency and also a damping of these oscillations. As in previous cases, the use of NCSTs strongly affects the ZF relaxation with an increasing effect with parameters $\gamma _k$, $\gamma _t$, $\gamma _n$ and $\sigma$. The effect of the smoothing is very strong even for $\sigma =0.1$. However, at short times the effect is not that large, and, anyway, much smaller than that of source terms when $\gamma _k$, $\gamma _t$ or $\gamma _n$ are above the ZF frequency.

4.3 Influence of NCSTs on nonlinear simulations

Once the effect of the NCSTs on the linear properties of unstable modes and ZFs has been studied we turn to study their effect and effectiveness in controlling numerical noise in nonlinear simulations in a stellarator configuration. We will use the same standard configuration of W7-X and study nonlinear simulations using the same ideal profiles as in § 4.2.1, with $T_{*} = 11.5\ \textrm {keV}$. This unrealistic temperature value is used just to make the simulations less expensive in computational resources, which allows us to run a set of them with different parameters at an affordable computational cost. We run simulations with adiabatic electrons. For these conditions, the growth rate of the most unstable mode, obtained from a reference linear simulation, is $\gamma _{\max }=8.7\times 10^{-4}\varOmega _*$, and the ZF oscillation frequency at half-radius is estimated as $\varOmega _{\textrm {ZF}}\approx 1.2\times 10^{-5}\varOmega _*$. This frequency is not calculated in a simulation but estimated using that obtained in § 4.2.2 and the known scaling of this frequency with temperature (Monreal et al. Reference Monreal, Sánchez, Calvo, Bustos, Parra, Mishchenko, Könies and Kleiber2017; Sánchez et al. Reference Sánchez, Calvo, Velasco, Medina, Alonso, Monreal and Kleiber2018). A square Fourier filter is used in all these simulations in order to reduce the numerical noise. The squared filter is defined by its width in poloidal and toroidal directions of the Fourier spectrum. Only modes with $-m_f < m < m_f$ and $-n_f < n < n_f$, with $n_f=150$ and $m_f=170$, are kept in the simulation. In addition to the square filter, a diagonal filter, aligned with the field line, is superimposed, which suppresses all modes with $|m - n /{\raise.3pt-\kern-5pt\iota} | > {\rm \Delta} m$, with ${\rm \Delta} m = 15$, because these modes have a parallel component which is too large. Note that only a period of the machine is simulated and, consequently, only modes with toroidal mode number being a multiple of the periodicity (five, in this case) are resolved.

In this study, we will pay attention to several quantities that allow one to measure the quality of the simulation: the signal-to-noise ratio, the volume-integrated zonal component and the non-zonal signal-to-signal ratio, which are defined as

(4.5)\begin{gather} \textrm{S2N} = \frac{\textrm{signal}}{\textrm{noise}} = \frac{\sum_s^{n_s} \sum_{m,n \in F} |\phi_{m,n}(s)|^2}{\sum_s^{n_s} \sum_{m,n \notin F} |\phi_{m,n}(s)|^2}, \end{gather}

(4.6)\begin{gather}\textrm{ZC} = \sum _s |\phi_{00}(s)|, \end{gather}

(4.7)\begin{gather}\textrm{NZ2S} = \frac{\text{non-zonal signal}}{\text{signal}} = \frac{\sum_s^{n_s} \sum_{m,n \in F, n\neq0} |\phi_{m,n}(s)|^2}{\sum_s^{n_s} \sum_{m,n \in F} |\phi_{m,n}(s)|^2}, \end{gather}

where $F$ represents the filter (intersection of the square and diagonal filters used in the simulation), $\phi _{\textrm {mn}}(s)$ is the Fourier component of the potential with poloidal and toroidal mode numbers $m$ and $n$, respectively (in $\theta ^*, \varphi$ coordinates, where $\varphi$ is the toroidal angle and $\theta ^*$ is the straight-field-line poloidal coordinate), which is summed in a finite number of flux surfaces, ${n_s}$, at which the potential spectrum is evaluated. The S2N ratio is computed assuming that all the power inside the filter is physically meaningful, while the power outside the filter corresponds to non-physical modes and is due to numerical noise.

In addition, we will look at the volume-integrated heat flux as a relevant physical quantity to study. In figure 13 the time evolution of the signal-to-noise ratio is shown for a set of simulations in which the NCSTs are used separately with different parameters $\gamma _k$, $\gamma _t$, $\gamma _n$ and $\sigma _{\textrm {QT}}, N_{\textrm {QT}}, f_{\textrm {QT}}$. The reference case, without any NCST, is shown in black colour in figures 13 and 14. It is clear from figure 13 that the heating source term has a small effect on the $S2N$ ratio for any value of $\gamma _t$, $\gamma _n$, while both Krook and weight smoothing have a significant effect. This means that a simulation in which only a source term is used, without including any noise control, will end up with the noise dominating the signal. The effect of both Krook and QT is an improvement of the quality of the simulation as measured by the S2N ratio, and this improvement increases with the size of parameters $\gamma _k$ and $\sigma _{\textrm {QT}}$. The Krook operator seems to be more efficient in increasing the S2N ratio. It is interesting to note that no convergence is found when these parameters are increased.

Figure 13.

Influence of the Krook operator (a), heating source (b) and weight smoothing for $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$ (c) on the signal-to-noise ratio (defined in (4.5)–(4.7)) in a W7-X equilibrium for different values of parameters.

Figure 14.

Influence of the Krook operator (a), heating source (b) and weight smoothing (for $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$) (c) on the zonal component in a W7-X equilibrium for different values of parameters. A volume-integrated ZC, as defined in (4.6), is used.

Another relevant quantity to look at is the zonal component (ZC). In figure 14 the volume-integrated ZC of the signal is shown versus time for several simulations with different parameters $\gamma _k$, $\gamma _t$, $\gamma _n$ and $\sigma _{\textrm {QT}}$. The reference case, without any NCST, is shown in black colour in all plots. It is clear that for noise-uncontrolled simulations, the ZC grows continuously in time without saturation. This feature was observed in previous nonlinear simulations with EUTERPE. This is related to the fact that the numerical noise directly contributes to the ZC. As for the S2N ratio, the heating sources have almost no effect on the ZC independently of the value of $\gamma _t$, $\gamma _n$. However, both Krook and weight smoothing have a significant effect on ZC, which increases as $\gamma _k$ and $\sigma _{\textrm {QT}}$ are increased, thus contributing to the stabilization of the ZC for long times. The effect of both tools is stronger over the ZC than over the S2N ratio and no clear saturation is observed in the range of parameters studied.

Now we look to the ratio of non-zonal component of the signal to the full signal, defined in (4.7). This quantity is shown versus time in figure 15 for the same set of simulations used in previous figures. In all cases, using a Krook operator, weight smoothing and heating sources, the NZ2S ratio increases as the parameters $\gamma _k$, $\gamma _t$, $\gamma _n$ and $\sigma _{\textrm {QT}}$ are increased. As in figures 13 and 14, the reference case, without any NCST, is shown in black colour. It is interesting to note in this figure that when no NCSTs are used, the NZ2S ratio saturates for times longer than $1.25 \times 10^{5} \varOmega ^{-1}$. We interpret this as a consequence of the simultaneous increase of noise and ZC for long times when no tools are used. This proves that a significant part of the ZC of the signal can be driven by the noise instead of being related to the nonlinear interaction of physical modes. When NCSTs are used, the NZ2S ratio increases up to a saturation higher level, which is close to the level just after the nonlinear saturation (at $t\sim 0.25$).

Figure 15.

Influence of the Krook operator (a), heating source (b) and weight smoothing for $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$ (c) on the ratio of the non-zonal component to the signal in a W7-X equilibrium for different values of parameters. The ratio NZ2S is defined in (4.7).

Finally, we turn to look at the heat flux, as a physically relevant quantity. Figure 16 shows the time trace of the volume-integrated heat flux, as defined in (3.2), for the same set of simulations studied before. As in previous figures, the reference case, without any NCST, is shown in black. It is clear that, if no NCSTs are used, the heat flux continuously decays after the initial overshot and the nonlinear saturation. This decay is the result of several effects. First, the instability has an effect on the profiles, which tends to relax the temperature gradient which drives the instability (ITG). On the other hand, the increase in the relative contribution of the uncorrelated numerical noise also has an effect on reducing the heat flux, which is produced by correlated fluctuations of density and velocity.

Figure 16.

Influence of the Krook operator (a), heating source (b) and weight smoothing for $N_{\textrm {QT}}=10, f_{\textrm {QT}}=\varOmega _*$ (c) on the volume-integrated heat flux in a W7-X equilibrium for different values of parameters.

When a Krook operator is included, the heat flux is stabilized to a value that increases with $\gamma _k$ and shows signs of some saturation for large enough values $\gamma _k>5\times 10^{-5}\varOmega _*$. Note that these are values close to $\gamma _{\max } /10$ as shown in § 4.1. As discussed there, for values $\gamma _k > \gamma _{\max } /10$ a significant effect of the Krook operator on the linear growth rate of unstable modes was observed in § 4.2.1, and the effect on the linear evolution of ZFs is strong even for smaller values of $\gamma _k$. Concerning the density and heating sources, we show in the figure the effect of switching on both density and heating sources separately. We can see that the heating is more effective in modifying the heat flux. This can be expected as it is the temperature gradient that drives the underlying instability (ITG). Then, restoring the ion temperature profile has a large effect on the heat flux. The effect of the density source on the heat flux is very small, and it is small in any case except that $\gamma _n$ is large as compared with the growth rate of unstable modes. It is observed that for $\gamma _t \ll \gamma _{\max } /10$ the heat flux slowly relaxes to a very small value. For values approaching $\gamma _t \sim \gamma _{\max } /10$ the slow decay of the heat flux after the initial overshoot is suppressed and the heat flux is stabilized, as shown by cases ‘$\gamma _t =5\times 10^{-5}$’ and ‘$\gamma _t =1\times 10^{-4}$’ in figure 16. Larger values can increase the heat flux level, it even increasing with time. Then, we can extract the rule of thumb that $\gamma _t$ should be $\gamma _t \leq \gamma _{\max } /10$ to stabilize the heat flux. Finally, the effect of weight smoothing on the heat flux is larger than that of density sources but smaller than that of Krook operator. A clear stabilizing effect is observed on the heat flux signal for values $\sigma _{\textrm {QT}} > 0.2$ with no significant effect on the heat flux level, which is a kind of saturation. If $\sigma$ is further increased, large heat flux spikes are observed without a significant increase in the average level (see the case $\sigma = 0.4$ in the figure). The saturation values of the heat flux are significantly smaller for the case of weight smoothing than for the case of using Krook operator.

After the comparison of several tools and their effect on the quality and physical results of the simulation, a short summarizing discussion is in order. As previously shown, the heating sources have a small effect on the simulation quality. Then, we can not rely on them for reducing the numerical noise. However, they can be used to stabilize the heat flux, and derived physical quantities, after the nonlinear saturation. It should be taken into account that the parameters $\gamma _t$ and $\gamma _n$ should be smaller than the relevant growth rate but, on the other hand, sufficiently large to sustain the equilibrium profiles without affecting the heat flux significantly. The Krook operator appears as an interesting tool for both the noise control and the heat flux stabilization. However, the strong dependency of physical quantities (growth rate and heat flux) with the $\gamma _k$ used makes it less adequate, particularly in the case of collisional simulations in which a Krook operator with values of $\gamma _k$ large enough as to contribute to the noise control can disturb significantly the effect of collisions. Another important drawback of using a Krook operator is the effect on the linear evolution of ZFs, which becomes important even at short times when $\gamma _k$ is large enough as to produce beneficial effects on the S2N ratio, as shown in figure 10. An energy and momentum conserving scheme as that implemented in ORRB5 could mitigate this negative feature, but implementing this for stellarator geometry is not as clear as in tokamak. In relation to the weight smoothing, it appears as an efficient tool for controlling the numerical noise and the secular zonal component which is driven by it. In addition, and similar to the Krook operator, weight smoothing has a stabilizing effect on the heat flux. As a positive point, we can remark that it has a very small effect on the linear ZF evolution at short times, which is expected to be the most relevant part for the nonlinear saturation of linearly unstable modes. As a drawback of this tool we can remark the lack of energy and momentum conservation, as in the process of smoothing the markers’ weights, both quantities are not conserved. The energy and momentum conserving implementation, proposed in (Donnel et al. Reference Donnel, Brunner, Gheller, Lanti, Ohana and Villard2019), has been tested in EUTERPE and the result is not very promising. When energy conservation is activated, then the effect of the weight smoothing is reduced up to almost not affecting the simulation, and its positive effect over the noise is lost. A partially conserving scheme is now under consideration. After these tests of the three NCS tools, it appears that for realistic cases, including collisions, a combination of heating sources with appropriate (minimum) parameters to sustain the profiles and a soft weight smoothing can be more adequate rather than using a Krook operator.

Before going to a more practical case, we next study how the heat flux value is affected when the radial limits of the physical domain of simulation are changed. This is important from a practical point of view, as reducing the radial extension of the computational domain can contribute to reducing the computational resources required for nonlinear global simulations which, in general, are very expensive.

4.4 Radially limited simulations

We have run a set of simulations with the same equilibrium from W7-X and profiles previously used in which the limits of the radial domain are changed. In all the simulations the characteristic scale lengths for density and ion temperature are $a/L_{\textrm {Ti}}=3$ and $a/L_n=1$, the simulations are carried out with adiabatic electrons, using the long-wavelength approximation, and a Krook operator, with $\gamma _k = 5\times 10^{-5} \varOmega _*$, is used to stabilize the simulations. The volume-integrated heat flux is shown versus time for this set of simulations in figure 17. The figure clearly shows that the heat flux level is not very sensitive to the location of the radial boundaries as far as the region of most unstable modes, which in this case are concentrated around half-radius, is captured in the simulation (solid time traces in the figure). However, the heat flux level is significantly reduced when the inner radial limit is increased above $r/a=0.4$, thus entering the location of unstable modes, and affecting them (dashed lines in the figure). In these simulations the instability parameter, $\eta _i=L_n/L_{\textrm {Ti}}$, is larger than 1 in the radial region $0.37 < r/a< 0.63$.

Figure 17.

Influence of the radial limits of the computational domain on the volume-integrated heat flux in a W7-X equilibrium for different locations of the radial boundaries. In all cases a Krook operator is used with $\gamma _k = 5\times 10^{-5} \varOmega _*$.

In practical cases, the criterion for limiting the radial domain may not be as clear as in this ideal case, as many different modes with different growth rates can coexist in different radial regions of the plasma, largely depending on the density and temperature profiles. However, this exercise gives us confidence on the limited effect of the special boundary region on the heat flux produced by modes located inside the simulation volume.

5 Application to a realistic W7-X case

Now we turn to study nonlinear stellarator simulations in a more realistic scenario. The purpose of this section is to demonstrate the use of the NCSTs in a stellarator configuration using relevant experimental plasma parameters of density and temperature and including collisions. For this purpose, we have selected an experimental program carried out in the standard magnetic configuration of W7-X (configuration Ref_170_EIM). We choose the W7-X experimental program 20181016_37, which has been extensively studied in (Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcuson, Alonso, Baldzuhn, Brandt, Brunner and Damm2020), is well diagnosed, and was used for power balance analysis. In this program, ECRH heating in the X2-mode was used, reaching an injected power of 5 MW. Much attention has been given in this program because cryogenic pellets were injected, thus increasing both the central density above $8 \times 10^{19}\ \textrm {m}^{-3}$ and peak ion temperature of $3\ \textrm {keV}$, and reaching record diamagnetic stored energy larger than $1.1\ \textrm {MJ}$. The focus of previous gyrokinetic analysis (Xanthopoulos et al. Reference Xanthopoulos, Bozhenkov, Smith, Plunk, Helander, Beidler, Alcuison, Grulke, Stechov and Alonso2020) has mainly been put on the external region of the plasma, $r/a>0.5$, where significant changes in profiles are produced in the post-pellet phase which lead to increased confinement (Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcuson, Alonso, Baldzuhn, Brandt, Brunner and Damm2020). However, we are now interested in testing the use of the NCSTs using experimentally realistic parameters with affordable resources. Computational resources required can be huge for the ion temperatures typical from the plasma edge in W7-X, because the wavenumber of the unstable modes increases as the temperature decreases and then, the spatial resolution has to be increased consequently. In addition, a wide spectrum of unstable modes (in terms of normalized wavenumber, $k_\bot \rho$) has been found in previous linear simulations in the outer region of W7-X (to be published), which further increases the computational resources required for simulating the plasma edge. This is the reason why we focus here in the innermost part of the plasma ($r/a<0.5$), where the ion temperature is above $1\ \textrm {keV}$. In this radial region of the plasma, during the post-pellet phase (seconds 1–2), the density profile is peaked and the ion temperature profile is quite flat, thus showing no ITG instability, which contributes to reducing turbulent transport in the core (Stechow et al. Reference Stechow, Grulke, Wegner, Proll, Alcuson, Smith, Xanthopoulos, Beidler, Beurskens and Bozhenkov2020). The stability of these profiles to ITG modes was confirmed with linear simulations (not shown here). In a later time after the pellet injection, the plasma relaxes to a poorer confinement state, with a flat density profile (with $n< 3.5 \times 10^{19}\ \textrm {m}^{-3}$) and peaked ion and electron temperature profiles at the same time (with central values $T_i \sim 1.5\ \textrm {keV}$ and $T_e \sim 3\ \textrm {keV}$), characteristic of ECH heated discharges in OP1.2, which features a situation with ITG instability (see Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcuson, Alonso, Baldzuhn, Brandt, Brunner and Damm2020). The profiles used in our simulations are shown in figure 18 and are characteristic of the normal confinement phase in that discharge,around t = 4s. We use analytical fits of the form $X = b[1-({r}/{a})^c]^d+e$, where $X$ represents $T_e$, $T_i$, or $n$. To investigate the role of collisions we run a simulation with a doubled density profile (see figure 18).

Figure 18.

Temperature and density profiles from W7-X program 20181016_37 (Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcuson, Alonso, Baldzuhn, Brandt, Brunner and Damm2020), for times around $t=4\ \textrm {s}$. The experimental density profile and another one with double the density are shown.

We will apply the NCSTs previously presented in a set of nonlinear simulations using these settings. We first run a linear simulation with these profiles without using any form of noise control. In this simulation the maximum linear growth rate is $\gamma _{\max }=4 \times 10^{-4} \varOmega _*$. Then, we have studied six different cases. First, as a reference case, we run a nonlinear simulation without accounting for collisions and without using any NCST. Next, we run another one including collisions, with a pitch angle scattering collision operator (Kauffmann et al. Reference Kauffmann, Kleiber, Hatzky and Borchardt2010; García-Regaña et al. Reference García-Regaña, Kleiber, Beidler, Turkin, Maaßberg and Helander2013). Then, a simulation including collisions and using a density profile with doubled density is considered in order to artificially enhance the influence of the collisions on the stabilization of heat flux and on the numerical noise, if any. The rest of the simulations all include collisions with the experimental density profile from figure 18 and using several NCSTs. The first one includes a Krook operator with moderate parameter $\gamma _k=5\times 10^{-6}\varOmega _*$ and also a heating source with $\gamma _h=5\times 10^{-6}\varOmega _*$. Next one is a simulation without Krook operator and using weight smoothing and heating sources, with $f_{\textrm {QT}}=1\varOmega _*, \sigma _{\textrm {QT}}=0.2, N_{\textrm {QT}}=10$ and $\gamma _h=5\times 10^{-6}\varOmega _*$. Finally, a simulation with weight smoothing and heating, using a less aggressive smoothing is considered, $f_{\textrm {QT}}=1\varOmega _*, \sigma _{\textrm {QT}}=0.1$, $N_{\textrm {QT}}=10$ and $\gamma _h=5\times 10^{-6}\varOmega _*$. The results are compared in figure 19, which shows the electrostatic energy, the S2N and NZ2S ratios and the volume-integrated heat flux for all these cases. First, we can highlight the weak influence of the collisions in all the measures studied. The results for cases including collisions (labelled with WCols), even with doubled density (case labelled with WColsNx2 in figure 19), are almost the same as without including them. In particular, the S2N and NZ2S ratios are very similar with or without collisions, but also the heat fluxes are almost the same. The S2N and NZ2S significantly improve both with Krook or weight smoothing. The influence of heating alone on S2N and NZ2S is very small (not shown), as expected, but it allows stabilizing the heat flux level. The heat flux level changes only slightly when weight smoothing is added, but the S2N quality improves significantly. Not much difference in the heat flux level is found between smoothing with $\sigma _{\textrm {QT}}=0.1$ or $\sigma _{\textrm {QT}}=0.2$. The Krook operator also improves the S2N ratio with a larger influence on the heat flux, which can change depending on the value of $\gamma _k$.

Figure 19.

Electrostatic field energy (a), S2N ratio (b), NZ2S ratio (c) and volume-integrated heat flux (d) for a set of simulations using W7-X standard configuration and experimental-like density and temperature profiles.

After the detailed study discussed in previous sections and the results presented in this one, we can conclude that stable reliable simulations with realistic parameters can be conducted using only heating sources to sustain the profiles and a soft weight smoothing in order to improve the simulation quality and without using a Krook operator. A Krook operator can also improve the S2N ratio; however, it has the drawback that can artificially affect collisions and, as demonstrated in previous sections, also affects the ZF dynamics, even in the short time. In addition, no clear saturation of the heat flux level is observed when $\gamma _k$ is increased. A heating source with a small $\gamma _h$ (a factor 100 below the maximum growth rate) parameter does not affect significantly neither the linear growth rate or the linear evolution of the ZF, as shown in § 4.2.2. Adding a soft weight smoothing, with $\sigma \sim 0.1$, contributes to improving the S2N ratio, stabilizing the ZC and the heat flux, while having a small effect on the short-time linear evolution of the ZFs. A good and steady S2N ratio, stable ZF component and heat flux can be reached with only these two elements.

6 Summary and conclusions

In this work, a comparison of global adiabatic-electron simulations carried out with the PIC codes ORB5 and EUTERPE in a tokamak configuration has been conducted. Good agreement was found in the linear properties of unstable modes between both codes.

Both growth rates and real frequencies were studied for different density and temperature profiles and a good agreement on the calculations of these quantities with both codes has been found. The EUTERPE has been benchmarked against ORB5 also in a nonlinear setting by means of simulations with adiabatic electrons in the same tokamak (CBC) equilibrium. The simulations have been conducted with the same numerical settings in both codes and, despite the differences between codes, a good agreement is found in the volume-integrated heat flux and heat conductivity provided that sufficient statistics is used.

In addition, a detailed study of several tools recently implemented in EUTERPE for the stabilization of simulations and the control of numerical noise has been carried out. Three tools have been thoroughly studied: heating sources implemented to sustain the density and temperature profiles, a Krook-type source term and also a smoothing of the markers’ weights which are aimed at reducing the numerical noise inherent to the PIC method and stabilizing the simulations. These tools have been tested in a tokamak (CBC) equilibrium and also in the standard configuration of the W7-X stellarator. The influence of using these tools on both linear (linear growth rate of unstable modes, linear relaxation of ZFs) and nonlinear (turbulent heat flux) properties has been studied in detail. The influence on the linear growth rate of unstable modes has been studied in tokamak and stellarator (W7-X) configurations while for the linear zonal flow evolution the focus is put on the stellarator and the specific features (low-frequency oscillation) appearing on it.

It is found that, as expected, the heating sources allow the sustaining of the profiles, with a limited effect on the linear properties and signal quality. Both Krook operator and smoothing affect the linear growth rate of modes and the linear relaxation of zonal flows. An increasing reduction of the linear growth rate is observed as the frequency of the Krook source term is increased, with moderate changes for values of $\gamma _k$ one order of magnitude below the linear growth rate. The weight smoothing has also an increasing effect on the linear growth rates which can be kept small for small values of the sigma parameter. While using a Krook source term produces a damping of ZFs at all time scales which is strong even for very small values of the frequency (well below the linear growth rate), the smoothing has a very limited influence on the time evolution of ZF for times shorter than time scales typical from the turbulence (turbulence decorrelation time or eddy turn-over time) in the order of tens of microseconds. For longer times the influence is also strong, however. Then, we can conclude that weight smoothing affects the ZF evolution less than the Krook operator at short times, which can be argued to be the most relevant time scales for turbulence saturation.

The influence of these tools in nonlinear simulations has also been studied. Several measures have been used to quantify the improvement of simulation quality. First, a signal-to-noise ratio has been used. In addition, the ratio of the non-zonal component of the signal to the total signal has been used. Special attention has been paid to the volume-integrated amplitude of the ZC of potential, which can be largely affected by the numerical noise. Finally, the influence of these tools on the volume-integrated heat flux has been studied. Using heating sources allows the sustaining of the density and temperature profiles and then helps in obtaining a stable heat flux signal, while its effect on simulation quality, as measured by S2N or ZF2S, is very small. Using either a Krook source term or the weight smoothing allows for improving of the simulation quality as measured by the S2N and ZF2S. Weight smoothing affects slightly the heat flux, while the influence of the Krook operator is stronger. The strong influence of Krook operator on the linear ZF damping at short times and the strong effect that it can produce on the turbulent heat flux makes the weight smoothing a preferred means for the control of numerical noise, particularly in collisional simulations in which the Krook source term can interfere and distort the effect of a more realistic collision operator. Stable simulations with improved (reduced) noise can be obtained using heating sources, with appropriate parameters, to sustain the profiles and weight smoothing for improving the simulation quality. The small damping of zonal flows at short times allows a proper treatment of ZF response. The weight smoothing presently implemented lacks of energy and momentum conservation, however. Work is in progress to develop a weight smoothing with conservation properties.

The reliability of a radially restricted simulation has been studied in a set of simulations in which different radial domains of computation are used instead of simulating the full volume. This is important in order to reduce the computation resources for a nonlinear simulation, which can be huge if the full radial domain is considered in realistic experimental conditions.

As an application of the tools for noise control and stabilization, a realistic case of W7-X has been studied. Radially restricted collisional nonlinear simulations using adiabatic electrons and realistic density and temperature profiles from a recent experimental program in W7-X have been carried out. The effect of NCSTs on the heat flux and quality measures has been studied in this practical case showing similar results to those observed in collisionless simulations with ideal density and temperature profiles carried out in the full radial domain. The turbulent heat flux, in this case, can also be stabilized using heating sources and weight smoothing, reaching a quasi-steady-state with reduced numerical noise.

To the best of our knowledge, these are the first nonlinear simulations carried out in a stellarator with a global PIC code, using adiabatic electrons, realistic experimental plasmas parameters and reaching a turbulent saturated steady-state, which have been reported.