1. Introduction: impurity ions in magnetised plasmas

Fully ionised classical plasmas are composed of electrons and ions of at least one, but often several, ion species. Any species other than the predominant ion species is commonly designated as an ‘impurity’. The transport of impurities is a subject of considerable relevance in toroidal fusion plasmas (Angioni Reference Angioni2021). An important transport mechanism in magnetised plasmas is turbulent advection of charged particles by the E×B drift velocity ${\boldsymbol {v}}_{E} = (1/|{\boldsymbol {B}}|^2) {\boldsymbol {E}} \times {\boldsymbol {B}}$ in the presence of fluctuating electric fields $\boldsymbol {E}(\boldsymbol {x}, t)$ and a static magnetic field $\boldsymbol {B}(\boldsymbol {x})$. The underlying theory and simulation of drift turbulence and instabilities in magnetised plasmas are covered by Scott (Reference Scott2021a,Reference Scottb).

Secondary ions of charge $q_j = Z_j e$ and particle density $n_j$ can be regarded as a ‘trace impurity’ with respect to turbulent transport when the charge density $Z_j e n_j \ll e n_e$ is much smaller than that of the electrons (index $e$), and are thus only passively advected by an E×B velocity (Lesur et al. Reference Lesur, Djerrou, Lim, Gravier, Idouakass, Moritz, Medina, Reveille, Drouot and Cartier-Michaud2020). This implies that the evolution of the governing electric field is solely determined by the quasi-neutral dynamics of the main ions ($i$) and electrons ($e$), for example as a result of instabilities or turbulence. Non-trace impurity ions with densities comparable to those of the dominant species can, however, contribute actively to the evolution of all field quantities, which necessitates a self-consistent model description including all ion species on the same level.

The nonlinear dynamics of turbulent advection of impurities in fusion edge plasmas has for example been studied with drift-fluid models (Naulin Reference Naulin2005; Priego et al. Reference Priego, Garcia, Naulin and Rasmussen2005; Dubuit et al. Reference Dubuit, Garbet, Parisot, Guirlet and Bourdelle2007; Futatani et al. Reference Futatani, del-Castillo-Negrete, Garbet, Benkadda and Dubuit2012; Bufferand et al. Reference Bufferand, Tamain, Baschetti, Bucalossi, Ciraolo, Fedorczak, Ghendrih, Nespoli, Schwander and Serre2019; Poulsen et al. Reference Poulsen, Rasmussen, Wiesenberger and Naulin2020), delta-F gyro-fluid models (Scott Reference Scott2005; Kendl Reference Kendl2012, Reference Kendl2014) and of impurities in an electron–positron plasma with a full-F gyro-fluid model (Kendl Reference Kendl2018). The E×B advection of impurities by interchange driven blobs has also been studied on the most fundamental particle level with a particle-in-cell model (Hasegawa & Ishiguro Reference Hasegawa and Ishiguro2017). Several delta-F and also more recently full-F gyro-kinetic studies have been devoted to combined neoclassical and turbulent impurity transport, often by employing quasi-linear or flux-surface averaged modelling (cf. exemplarily Garcia-Regana et al. (Reference Garcia-Regana, Barnes, Calvo, Parra, Alcuson, Davies, Gonzales-Jerez, Mollen, Sanchez and Velasco2021), and the review by Angioni (Reference Angioni2021) and references therein).

In comparison with drift-fluid models, gyro-fluid models have several advantages when it comes to the modelling of impurities. As was recently pointed out by Poulsen et al. (Reference Poulsen, Rasmussen, Wiesenberger and Naulin2020), drift-fluid models suffer from a species asymmetry: two separately modelled, identical ion species are not equivalent to a single ion species of summed density. This problem is not present in gyro-fluid models, where all species are modelled by an equivalent set of equations, as we will explicitly show in § 2. Further, drift-fluid equations are strongly coupled in the time domain through the polarisation drift. This leads to an intricate inversion problem in numerical implementations (Poulsen et al. Reference Poulsen, Rasmussen, Wiesenberger and Naulin2020), in particular if an implicit time step is used. Bufferand et al. (Reference Bufferand, Tamain, Baschetti, Bucalossi, Ciraolo, Fedorczak, Ghendrih, Nespoli, Schwander and Serre2019) uses additional approximations to break this coupling. Gyro-fluid models are free from this problem altogether as the species coupling via the polarisation equation is free from time derivatives. The model is thus much easier to implement. The last advantage of gyro-fluid models lies in their inherent capability to treat finite Larmor radius (FLR) effects consistently to high order. Often, impurity ion species have large gyro-radii through their large mass and possible incomplete ionisation, which may break the assumption $\rho _i k_\perp \ll 1$, where $\rho _i = \sqrt {m_i T_i} / (Z_i e B)$ is the gyro-radius of ions with mass $m_i$ and Temperature $T_i$, and $k_\perp$ is the typical fluctuation wavenumber.

Delta-F gyro-fluid or gyro-kinetic models for impurity advection transport studies cannot consistently handle spatially localised trace or non-trace impurities (Kendl Reference Kendl2014), because these assume the presence of an ‘infinitely large’ reservoir of a stationary background component ${F_0}$, while dynamically treating only small fluctuations $\delta F$ of the particle distribution functions $F = {F_0} + \delta F$. A full-F formulation including the total impurity distribution function $F_j$ is at least necessary for the trace impurity component (while the main species could be still approximated by a delta-F formalism). Moreover, for non-trace impurities a consistent formulation is required that couples all species in full-F formalism through the quasi-neutrality constraint in the polarisation equation.

In the following, the applicability of a self-consistent full-F gyro-fluid model on non-trace impurity advection driven by interchange unstable blob modes is demonstrated. The transition from trace to non-trace cases, involving consistent coupling to the evolution of all field quantities, is quantified.

2. Full-F multi-species gyro-fluid model of interchange blob transport

Full-F gyro-kinetic models evolve the full five-dimensional gyro-kinetic distribution functions $F_s(\boldsymbol {x}, v_{\perp }, v_{\parallel }, t)$ for all charged particle species, without referring to the computationally greatly simplifying limit of treating only small fluctuations $\delta F (\boldsymbol {x}, v_{\perp }, v_{\parallel }, t)$ on a stationary equilibrium background ${F_{s0}} = F_s -\delta F_s$. Gyro-fluid models can be obtained by evolving velocity moments of $F_s$ with an appropriate closure. Full-F gyro-fluid models have been derived first by construction from a postulated gyro-fluid Lagrangian (Strintzi & Scott Reference Strintzi and Scott2004, Reference Strintzi and Scott2005) and later by a moment-based derivation from full-F gyro-kinetics (Madsen Reference Madsen2013). Full-F isothermal gyro-fluid simulations based on the model by Madsen (Reference Madsen2013) were first applied to the case of finite ion temperature, high amplitude plasma blobs driven by the interchange instability in an inhomogeneous magnetic field in two dimensions by Wiesenberger, Madsen & Kendl (Reference Wiesenberger, Madsen and Kendl2014) and also three dimensions by Kendl (Reference Kendl2015). Temperature dynamics, through higher-order moments and variable FLR effects, was introduced by Held et al. (Reference Held, Wiesenberger, Madsen and Kendl2016).

Subsequently, also drift wave-interchange turbulence, zonal flow and vortex dynamics were addressed with this full-F gyro-fluid model (Held Reference Held2017; Held et al. Reference Held, Wiesenberger, Kube and Kendl2018; Held, Wiesenberger & Kendl Reference Held, Wiesenberger and Kendl2019). From these studies it has emerged that in the assumed long-wavelength (‘low-k’) approximation $\rho _i k_{\perp } \ll 1$ of Strintzi & Scott (Reference Strintzi and Scott2004, Reference Strintzi and Scott2005) and Madsen (Reference Madsen2013), is, however, missing a physically relevant FLR effect in the polarisation density that appears in the quasi-neutrality constraint (Held, Wiesenberger & Kendl Reference Held, Wiesenberger and Kendl2020). Arbitrary wavelength (‘full-k’) polarisation, which consistently incorporates the additional FLR physics via exact and Padé-approximated closure relations, was introduced to full-F gyro-fluid theory by Held et al. (Reference Held, Wiesenberger and Kendl2020) and numerical studied in Held & Wiesenberger (Reference Held and Wiesenberger2022). The latter numerically study shows that ‘full-F full-k’ effects result in highly coherent, long-lived and decelerating blob structures. The most general ‘full-F full-k’ form of the polarisation equation is, however, computationally more demanding: Padé-approximated numerical solvers based on Held et al. (Reference Held, Wiesenberger and Kendl2020) can be more readily constructed for isothermal gyro-fluid models with a single ion species (Held & Wiesenberger Reference Held and Wiesenberger2022) while requiring only a slight computational overhead in comparison with low-k or Delta-F simulations. The generalisation to multiple ion species will be addressed in future works.

The present investigation of impurity advection is therefore based on a conservative multi-species full-F gyro-fluid model in the long-wavelength (low-k) approximation, which is still a unique approach that has so far not been studied.

Tokamak plasma edge dynamics regulates the generation and distribution of impurities (which basically are all non-fuel particles) in a fusion plasma, but is itself affected by the presence of these ion species. At the plasma core impurity accumulation degrades fusion performance, but at the edge impurity effects may mitigate catastrophic heat load to plasma facing components. In L-mode operation, experimental evidence attributes up to half of the heat and particle transport from the scrape-off layer (SOL) towards the first wall to ‘blob’ structures: dense and hot plasma filaments, elongated along the magnetic field, which detach from the plasma and propagate radially outwards (Krasheninnikov, Smolyakov & Kukushkin Reference Krasheninnikov, Smolyakov and Kukushkin2020).

Perturbation amplitudes – in particle densities and temperatures – of such structures, compared with the background quantities in the SOL of fusion experiments, can be well above unity. This defies any delta-F type models based on the separation of background and fluctuating quantities. Consistent blob evolution and transport has previously mostly been studied with global fluid models and codes, such as STORM/BOUT$++$ (Easy et al. Reference Easy, Militello, Omotani, Dudson, Havlickova, Tamain, Naulin and Nielsen2014; Hoare et al. Reference Hoare, Militello, Omotani, Riva, Newton, Nicholas, Ryan and Walkden2019), GBS (Ricci et al. Reference Ricci, Halpern, Jolliet, Loizu, Mosetto, Fasoli, Furno and Theiler2012), HESEL (Nielsen et al. Reference Nielsen2016; Poulsen et al. Reference Poulsen, Rasmussen, Wiesenberger and Naulin2020), TOKAM3X (Tamain et al. Reference Tamain, Bufferand, Ciraolom, Colin, Galassi, Ghendrih, Svhwander and Serre2016), SOLT (Russel et al. Reference Russel, Myra, D'Ippolito, LaBombard, Hughes, Terry and Zweben2016) or GRILLIX (Ross et al. Reference Ross, Stegmeir, Manz, Groselj, Zholobenko, Coster and Jenko2019), and the importance of FLR effects has been demonstrated with a delta-F gyro-fluid model (Madsen et al. Reference Madsen, Garcia, Larsen, Naulin, Nielsen and Rasmussen2011). Recently, a full-F gyro-kinetic code (Myra et al. Reference Myra, Ku, Russell, Cheng, Keramidas-Charidakos, Parker, Churchill and Chang2020) and a full-particle particle-in-cell (PIC) code (Hasegawa & Ishiguro Reference Hasegawa and Ishiguro2017) have been applied to the analysis of blob filament propagation, the latter also including a secondary ion species.

In the following, the full-F gyro-fluid model based on Madsen (Reference Madsen2013) is used, which is computationally much more economical compared with a gyro-kinetic or N-body particle description, but still retains important FLR effects, which are not present in fluid models. The main novel focus here, however, is on the inclusion of an additional (impurity) ion species, which is readily incorporated within a full-F gyro-fluid (or gyro-kinetic) model.

The radial blob transport properties are determined by two-dimensional (2-D) advection in the $(x,y)$ drift plane perpendicular to a magnetic field $\boldsymbol {B}$. Parallel evolution of blobs, which are initially localised along the magnetic field, can contribute significantly to the overall dynamics, but is here neglected in order to highlight the convective part of the multi-species dynamics. This physically corresponds to the case of strongly elongated filaments with small parallel gradients, which obey a quasi-2-D drift dynamical evolution. For the same reason of restriction to 2-D advective dynamics only, the otherwise important collisional interactions between the impurities and the main plasma ions and electrons are also neglected here.

The set of isothermal full-F gyro-fluid equations by Wiesenberger et al. (Reference Wiesenberger, Madsen and Kendl2014) is here extended to describe the evolution of gyro-centre densities $N_s$ of all charged particle species, with $s \in (e,i,j)$ specifically including electrons ($e$), main ions ($i$) and an impurity ion species ($j$), which are coupled by quasi-neutrality through the long-wavelength form of the gyro-fluid polarisation equation.

The resulting equations for the evolution of the gyro-centre densities $N_s(\boldsymbol {x},t)$ for each species are

(2.1)\begin{equation} \partial_t N_s + \boldsymbol{\nabla}\boldsymbol{\cdot} \left((\boldsymbol{U}_{\boldsymbol{E},\boldsymbol{s}} + \boldsymbol{U}_{\boldsymbol{\nabla} \boldsymbol{B},\boldsymbol{s}})N_s\right) ={-}\nu \boldsymbol{\nabla}_{{\perp}}^4 N_s, \end{equation}

where $\boldsymbol {U}_{\boldsymbol {E},\boldsymbol {s}} = (\boldsymbol {e}_{\boldsymbol {z}} \times \boldsymbol {\nabla }\boldsymbol {\psi }_{\boldsymbol {s}})/B$ and $\boldsymbol {U}_{\boldsymbol {\nabla } \boldsymbol {B},\boldsymbol {s}} = T_s (\boldsymbol {e}_{\boldsymbol {z}} \times \boldsymbol {\nabla } \ln B)/(q_sB)$ with $T_s$ the constant species temperature and $q_s = Z_s e$ the species charge. Electrons are here assumed to have a vanishing gyro-radius compared with the ion species, so that the electron particle density $n_e = N_e$ coincides with the electron gyro-centre density, and all FLR effects through gyro-averaging operators on electrons are neglected.

The inhomogeneous static background magnetic field strength $B(x)= B_0 / (1+x/R)$ is assumed to only depend on the radial ($x$) direction, where $R$ would correspond to the major radius of a toroidal plasma. We then set $\boldsymbol {B} = \boldsymbol {B}(\boldsymbol {x}) \boldsymbol {e}_{\boldsymbol {z}}$, identifying $\boldsymbol {e}_{\boldsymbol {z}}$ as the magnetic unit vector. The radial gradient in this magnetic field provides the interchange drive of perturbations with a co-aligned local or global pressure gradient through the term $\boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {U}_{\boldsymbol {\nabla } \boldsymbol {B},\boldsymbol {s}} \boldsymbol {N}_{\boldsymbol {s}}) = - T_s \partial _y N_s / (q_s R B_0)$. For more details see Held et al. (Reference Held, Wiesenberger, Madsen and Kendl2016).

The hyperviscosity term with small parameter $\nu$ is added for numerical stability.

The electric plasma potential $\phi (\boldsymbol {x},t)$ appears in the full-F gyro-fluid potential

(2.2)\begin{equation} \psi_s = \varGamma_s \phi - \frac{1}{2} \frac{m_s}{q_s} \left( \frac{\boldsymbol{\nabla}_{{\perp}} \phi}{B} \right)^{2} \end{equation}

where the gyro-averaging operator $\varGamma _s = ( 1- \frac {1}{2} ({m_s T_s}/{q^2B^2})\boldsymbol {\nabla }_{\perp }^2 )^{-1}$ is expressed in Padé approximation, and the second term is the full-F correction to the potential through the E×B energy in its long-wavelength approximation.

Because of the small electron mass and vanishing electron gyro-radius ($\varGamma _e = 1$) in comparison with ions, $\varPsi _e = \phi$ can be assumed.

The electric potential $\phi$ is determined through the quasi-neutral nonlinear polarisation equation in full-F long-wavelength form

(2.3)\begin{equation} \sum_s \left[ \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \frac{m_s N_s}{B^2} \boldsymbol{\nabla}_{{\perp}} \phi \right) \right] ={-} \sum_s \left[ q_s \varGamma_s N_s \right] . \end{equation}

This multi-species form of the polarisation equation shows that any additional charged particle species is formally simply added as separate terms in the sums, where the gyro-centre density $N_s$ of each species is evolved by an own equation of (2.1).

Equation (2.3) is the expression of quasi-neutrality between particle species $\sum _s q_s n_s = 0$ where the particle density $n_s$ is transformed from the gyro-centre density $N_s$ via

(2.4)\begin{equation} n_s = \varGamma_s N_s + \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \frac{m_s N_s}{ q_s B^2} \boldsymbol{\nabla}_{{\perp}} \phi \right). \end{equation}

Note that neglecting electron mass $m_e = 0$ we obtain $n_e = N_e$, as mentioned above.

By using the abbreviation $\epsilon _s = m_s N_s / B^2$, the left-hand side of (2.3) carrying the polarisation densities can, in this full-F low-k form, be re-written as

(2.5)\begin{equation} \sum_s \left[ \boldsymbol{\nabla}\boldsymbol{\cdot} ( \epsilon_s \boldsymbol{\nabla}_{{\perp}} \phi ) \right] = \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \left[ \sum_s \epsilon_s \right] \boldsymbol{\nabla}_{{\perp}} \phi \right) \equiv \boldsymbol{\nabla}\boldsymbol{\cdot} ( {\bar{\epsilon}} \,\boldsymbol{\nabla} \phi ) = Q. \end{equation}

The polarisation equation (2.3) therefore again takes the same form of a generalised Poisson equation with a spatially dependent dielectric coefficient ${\bar {\epsilon }} (\boldsymbol {x},t)$ and source function $Q(\boldsymbol {x},t) = - \sum _s [ q_s \varGamma _s N_s ]$, and can be treated with the same numerical solvers as implemented for codes with a single ion species.

It is important to realise the inherent species symmetry in the gyro-fluid approach. The density equation (2.1) for separate ion species differs only by the species index $s$, that is, different mass, charge and temperature. The equations for two separate but identical ion species $N_{1s}$ and $N_{2s}$ can be summed to yield

(2.6)\begin{equation} \partial_t (N_{1s}+N_{2s}) + \boldsymbol{\nabla}\boldsymbol{\cdot} \left((\boldsymbol{U}_{\boldsymbol{E},\boldsymbol{s}} + \boldsymbol{U}_{\boldsymbol{\nabla}\boldsymbol{B},\boldsymbol{s}}) (\boldsymbol{N}_{\boldsymbol{1s}}+\boldsymbol{N}_{\boldsymbol{2s}})\right) ={-}\nu \boldsymbol{\nabla}_{{\perp}}^4 (N_{1s}+N_{2s}), \end{equation}

which is equivalent to the equation for a single ion species of density $N_s = N_{1s} + N_{2s}$. Further, note that the coupling between the species is given by the polarisation equation (2.3), which is free from a time derivative. This equation can thus be solved in each time step separately in a numerical implementation. As was pointed out in § 1 this separates gyro-fluid models from drift-fluid models.

The Bohm normalisation (cf. Wiesenberger et al. Reference Wiesenberger, Madsen and Kendl2014; Held et al. Reference Held, Wiesenberger, Madsen and Kendl2016) is used to renormalise the equations, based on the gyration frequency $\varOmega _i =e B_0 / m_i$ of the main ions with mass $m_i$, and the drift scale $\rho _0 = \sqrt {m_iT_e}/(eB_0)$ with respect to the electron temperature $T_e$, which are related through the sound speed $c_0 = \rho _0 \varOmega _i$. The densities are normalised individually by choosing values $N_{s0}$ while the potential is normalised with $e/T_e$.

The dimensionless parameters of the model are then given by

(2.7a–d)\begin{equation} \tau_s =\frac{T_s} { Z_s T_e} \quad \mu_s = \frac{m_s}{Z_s m_i} \quad a_s = \frac{Z_s N_{s0}}{N_{e0}} \quad \kappa = \frac{\rho_0}{R}.\end{equation}

The constant temperatures enter through the non-dimensional ratio $\tau _s$, so that always $\tau _e = -1$, while $\tau _i$ and $\tau _j$ are variables. Further, the normalised mass ratio $\mu _s$ is introduced as well as the density ratios $a_s$. The magnetic field gradient driven interchange terms appear with a ‘curvature’ parameter $\kappa$. The dimensionless set of equations we discuss in the following then reads

(2.8)\begin{gather} \partial_t N_s + \frac{1}{B} \left \{ \psi_s, N_s \right \} = N_s \kappa \partial_y \psi_s + \tau_s \kappa \partial_y N_s -\nu \boldsymbol{\nabla}_{{\perp}}^4 N_s. \end{gather}

(2.9)\begin{gather}\sum_s \left[ \boldsymbol{\nabla}\boldsymbol{\cdot} \left( \frac{a_s\mu_s N_s}{B^2} \boldsymbol{\nabla}_{{\perp}} \phi \right) \right] ={-} \sum_s \left[ a_s \varGamma_s N_s \right], \end{gather}

with

(2.10a,b)\begin{equation} \psi_s = \varGamma_s \phi - \frac{1}{2}\mu_s \left(\frac{\boldsymbol{\nabla}_\perp\phi}{B}\right)^2 \quad \varGamma_s = \left(1-\frac{1}{2}\tau_s \mu_s \boldsymbol{\nabla}^2_\perp\right)^{{-}1}. \end{equation}

For boundary conditions we choose Dirichlet in $x$, with $N_s = 1$, $\phi = \psi _s = 0$, and periodic boundaries in the $y$ direction. For the simulations presented in §§ 4 and 5 the computational domain is chosen large enough that the specific choice of boundary conditions does not influence the results.

3. Multi-species treatment in FELTOR code

Solutions of two- and multiple ion species formulations of (2.1) and (2.3) are calculated within the FELTOR (Full-F ELectromagnetic code in TORoidal geometry) C$++$ numerical framework (Wiesenberger & Held Reference Wiesenberger and Held2022). We choose a discontinuous Galerkin (dG) scheme for spatial discretisations and an adaptive explicit embedded Runge–Kutta method as time integrator. Beside common design goals – fast development, quick implementation of model equation sets and high speed in performing the necessary calculations – FELTOR especially aims at reproducibility (Wiesenberger et al. Reference Wiesenberger, Einkemmer, Held, Gutierrez-Milla, Saez and Iakymchuk2019), being a freely available, modular scientific software package. In particular, we make the C$++$ code and Python data-analysis scripts available that can be used to reproduce the simulations and plots in Wiesenberger (Reference Wiesenberger2022).

The performance of the present multi-species model implementation is given by three main factors: the runtime for (i) the right-hand side of (2.8) and the operations for the time stepper, (ii) the solution of the Helmholtz-type elliptic operators $\varGamma _s N_s$ and $\varGamma _s \phi$ and (iii) the solution of the elliptic equation (2.3) for $\phi$. Note here that the inversion of the $\varGamma _s$ only plays a role if $\tau _s\neq 0$, which in this paper is the case only in some of the simulations in § 5. Furthermore, two $\varGamma _s$ need to be inverted per ion species while the polarisation equation needs to be inverted only once per time step.

The implementation of the right-hand side of (2.8) and any single- or multi-step time stepper can be done with (almost) perfectly parallelisable matrix–vector multiplications and trivially parallel vector–vector operations. These operations are extremely fast and, per species, amount to approximately a single iteration of a conjugate gradient solver. If $\tau _s\neq 0$ the right-hand side is thus negligible performance-wise compared with the inversion of $\varGamma _s$. For $\tau _s =0$ only a high number of species ($10$ say) will have an impact on the resulting simulation time, while the solution of the polarisation equation remains the dominant factor.

The inversion of a single $\varGamma _s$ typically takes approximately $1/5$ of the time to solve the polarisation equation. For even two species with $\tau _s\neq 0$ the cumulative time to solve all $\varGamma _s$ can thus become equal to that of the polarisation equation. If the amount of impurity species is increased further, the inversion of $\varGamma _s$ may even become the dominant performance indicator.

Both the inversion of $\varGamma _s$ and the polarisation equation require an efficient elliptic solver and in FELTOR we use the same solver for both. Thus, both in the case $\tau _s = 0$ as well as $\tau _s\neq 0$ the elliptic solver is the main bottleneck for performance. In order to illustrate the performance of elliptic solvers we study the polarisation equation (2.5) in an isolated setting. The inversion of $\varGamma _s$ is analogous. We solve

(3.1)\begin{equation} -\boldsymbol{\nabla}\boldsymbol{\cdot}( \chi \boldsymbol{\nabla}_\perp \phi) = \rho \end{equation}

for $\phi$, given $\rho$ and $\chi$. We use (local) dG methods Cockburn & Shu (Reference Cockburn and Shu1998) to discretise $\partial _x \rightarrow D_x$, where $D_x$ is a block-sparse matrix. We then have

(3.2)\begin{gather} M \phi = \rho, \end{gather}

(3.3)\begin{gather}M = D_x^T \chi D_x + D_y^T \chi D_y + J . \end{gather}

We here see that $M$ is self-adjoint, which means that we can use a conjugate gradient (CG) solver. A CG solver is an example of a matrix-free solver, that is a solver for $M x = b$ that does not require access to the elements of the matrix $M_{ij}$. We use this property in FELTOR to implement $M$ in a stencil-like fashion such that $M$ fits entirely in the cache and thus memory is saved. Matrix-free solvers also have the advantage that the matrix $M$ does not need to be assembled in every time step. Thus an expensive matrix–matrix multiplication (as in $D_x^{T} \chi D_x$) in each time step can be avoided. Recall here that $\chi$ in (3.1) or equivalently ${\bar {\epsilon }}$ in (2.5) changes at each time step. The downside is that, since the matrix elements $M_{ij}$ are unknown, the preconditioner for $CG$ also needs to be matrix free. Our best guess is thus a diagonal preconditioner

(3.4)\begin{equation} M^{{-}1} \approx \frac{1}{\chi}. \end{equation}

We manufacture a solution to (3.1)

(3.5)\begin{gather} \chi = 1 + A \sin (x)\sin(y) \end{gather}

(3.6)\begin{gather}\rho = 2 \sin(x) \sin(y) (A\sin(x)\sin(y) +1)- A\sin^2(x)\cos^2(y) - A\cos^2(x)\sin^2(y) \end{gather}

(3.7)\begin{gather}\phi = \sin(x)\sin(y) \end{gather}

with $A = 0.9$ and solve on the domain $[0,{\rm \pi} ]\times [0,2{\rm \pi} ]$ for Dirichlet boundary conditions in $x$ and periodic in $y$. As resolution we choose $3$ polynomial coefficients and use $512$ cells in $x$ and $y$. Each vector thus consists of $1536\times 1536$ points or approximately $20$ MB of data. The initial guess is zero. This problem reflects to a high degree the condition and behaviour of the inversion problem that has to be solved in every time step for the blob simulations in §§ 4 and 5. This is true even though in the simulation an initial guess based on previous time step solutions is available.

We compare various solver methods in figure 1.

Figure 1. Performance of different solver methods in FELTOR for (3.1) on a NVIDIA GeForce RTX 3090. Note the logarithmic scale on the $y$-axis.

The solver with the longest time to solution is a LGMRES (loose generalised minimal residual method) solver with diagonal preconditioner (3.4). It takes almost 10 times as long as the preconditioned CG solver. This is due to both a much higher iteration number as well as more time for each iteration. The next faster solver is a biconjugate gradient stabilised (BICGSTABl) scheme that takes almost as long as an unpreconditioned CG solver. With these examples it becomes apparent that our symmetric/self-adjoint discretisation (3.3) is highly beneficial. One should not try to discretise a non-symmetric form of (3.1) like $-(\chi \varDelta + \boldsymbol {\nabla } \chi \boldsymbol {\cdot }\boldsymbol {\nabla }) \phi = \rho$ because non-symmetric solvers like LGMRES or BICGSTABl are (potentially a lot) slower than CG.

The highest impact on performance has, however, our nested iterations ‘full approximation scheme’ (FAS). This is a type of matrix-free multigrid scheme where on each grid a pre-conditioned conjugate gradient (PCG) solver is used to entirely solve the problem. First, an initial guess is projected down to the coarsest grid and the matrix problem is solved with PCG (FELTOR in fact allows any solver but PCG works best for symmetric problems, as shown above). The solution is then interpolated onto the next finer grid and used as an initial guess for PCG there. This process is then repeated until the finest grid is reached. We typically use 3 grids. As can be seen in figure 1, the performance improvement over a simple PCG is dramatic, with more than factor 10. We also note that the performance of our FAS scheme can be further improved by avoiding scalar products (global communication) on the coarsest grids. This can be done by evaluating the stopping criterion for PCG only every tenth iteration, for example. This gives another slight performance benefit in the optimised FAS setting.

Finally, we note that increasing the amplitude in the manufactured solution (3.5) increases the number of iterations of the solvers. For example, setting the amplitude to $A=0.99$ increases the runtime of our FAS scheme by approximately 10 %. This reflects a somewhat indirect influence of additional species on performance compared with only a single ion species. The change in amplitude corresponds to a change of $\bar {\epsilon }$ in the polarisation (2.5). As long as $a_j\mu _j \ll 1$, the impact is low but for higher concentrations as in § 5 this factor becomes relevant.

4. Background impurity effects on interchange blob evolution

The first test scenario of the multi-species full-F code assumes a constant background impurity density with $a_j < a_i$, and studies the impact of relative impurity density $a_j$ and mass $\mu _j$ on blob propagation. This corresponds to a transition from trace to non-trace impurity conditions. The initial Gaussian seeded perturbation is here assumed to have the same relative impurity concentration as the overall background plasma.

In order to focus on full-F impurity effects in the present model, the first simulations presented here are for cold main and impurity ions with $\tau _i=0$ and $\tau _j=0$. This excludes FLR effects and associated poloidally asymmetric propagation (Madsen et al. Reference Madsen, Garcia, Larsen, Naulin, Nielsen and Rasmussen2011), and facilitates a better comparability of the maximum radial blob velocity $v_{\max }^x$ as a measure for the impurity concentration effect on blob transport.

From linearisation of the gyro-fluid model (2.1) and (2.3), the growth rate $\gamma$ of the interchange unstable system can be obtained as a function of the initial Gaussian blob width $\sigma$ with relative amplitude $A$.

This results in an analytical scaling law for the maximum blob velocity via a scaling analysis of the vorticity equation as

(4.1)\begin{equation} v_{\max}^x = \gamma \sigma \sim \sqrt{ \frac{A \kappa \sigma (1+{\bar{\tau}})}{\bar{\mu}}}, \end{equation}

with an effective pressure ratio ${\bar {\tau }} = \sum a_s \tau _s$ and mass density ratio ${\bar {\mu }} = \sum _s a_s \mu _s$. For cold ions and impurities ${\bar {\tau }} = 1$, so that the impurity concentration effect only enters through the effective mass density as $v_{\max }^x \sim 1/\sqrt {\bar \mu }$.

We here point out that there exists a second ‘energy-limited’ (or ‘linear’) regime where the maximum blob velocity depends linearly on the amplitude (Kube, Garcia & Wiesenberger Reference Kube, Garcia and Wiesenberger2017; Wiesenberger et al. Reference Wiesenberger, Held, Kube and Garcia2017; Held & Wiesenberger Reference Held and Wiesenberger2022). This regime occurs when the blobs’ initial energy is insufficient to reach the maximum velocity of the square root scaling even if all of the initial energy is converted to kinetic energy. The blobs in the following study are outside of this regime.

In the nonlinear simulations, the instantaneous velocity $v^x(t)$ is measured through the change in position of the evolving blob's centre of mass over time. After a nearly linear acceleration phase, the blob velocity typically reaches a maximum, before it slows down again mostly due to nonlinear break up. The maximum velocity $v_{\max }^x$ can so be determined in simulations from the resulting computed function $v^x(t)$.

We define the $x$-component of the gyro-centre of mass position and the gyro-centre of mass of species $s$ as

(4.2)\begin{gather} X_{s} = \frac{1}{M_s} \int {\rm d}V \,x (N_s -1), \end{gather}

(4.3)\begin{gather}M_s = \int {\rm d}V (N_s-1). \end{gather}

It is straightforward to see that the gyro-centre mass $M_s$ is a time invariant $\dot {M} = 0$ neglecting dissipation $\nu = 0$ in (2.1). The time derivative of $X_s$ yields

(4.4)\begin{align} \frac{{\rm d}}{{\rm d} t} X_s & = \frac{1}{M_s} \int {\rm d}V \,x \dot{N}_s = \frac{1}{M} \int {\rm d}V \,N_s \partial_y \psi_s /B \end{align}

(4.5)\begin{align} & =\frac{1}{M} \int {\rm d}V \,N_s \partial_y ( \varGamma_s \phi - 0.5 \mu_s (\boldsymbol{\nabla} \phi)^2/B^2) /B \end{align}

(4.6)\begin{align} & =\frac{1}{M} \int {\rm d}V \,\varGamma_s N_s \partial_y \phi /B + \partial_y \phi /B \boldsymbol{\nabla}\boldsymbol{\cdot} \mu_s N_s (\boldsymbol{\nabla}_\perp \phi)/B^2) \end{align}

(4.7)\begin{align} & =\frac{1}{M_s} \int {\rm d}V \,n_s u_E^x, \end{align}

where $u_E^x = \partial _y \phi / B$ is the $x$-component of the E×B velocity and $n_s$ is the particle density of species $s$. We use (2.1) and (2.4) and neglect terms from surface integration. The gyro-centre of mass velocity is thus directly related to the radial particle transport of species $s$. For the following analysis we use the main ion centre of mass velocity as $v^x(t)$.

The maximum centre of blob mass velocity is obtained by a series of simulations as a function of variable relative impurity concentrations $a_j$ and impurity mass $\mu _j$. All other parameters are here fixed by assuming blob dimensions that are representative of structures typically encountered in the SOL of present medium sized tokamaks: initial relative perturbation amplitude $A = \Delta N_i / N_{i0} = 1$, and blob width $\sigma = 10$ in units of $\rho _0$. The curvature parameter $\kappa = 0.000457$ is chosen to represent a parameter set typical for the low-field side SOL of, for example, ASDEX Upgrade experiments. The main ions are assumed to be singly charged ($Z_i = 1$).

The impurity concentration is varied in the range $10^{-3} \leqslant a_j \leqslant 0.5$, and the relative mass within the range $2 \leqslant \mu _s \leqslant 20$. For example, a W$^{30+}$ tungsten impurity with a concentration of $10^{-4}$ relative to deuterium main ions would correspond to parameters $a_w = Z_w N_{w0} / N_{e0} = 3 \times 10^{-3}$ and $\mu _w = m_w / (Z_w m_i) = 3$, and a N$^{2+}$ nitrogen seeded SOL impurity with concentration $0.05$ would correspond to $a_N = 0.1$ and $\mu _N = 7$. A value of $\mu _j=20$ could for example correspond to singly charged argon ions. The present model could be applied with an arbitrary number of ion species (with accordingly raised computational costs), which could also mean different ionisation states of a specific element. In the following, only one impurity species is included.

The nonlinear simulation results for $v_{\max }^x$ as functions of $a_j$ and $\mu _j$ are shown in figure 2 as symbols, together with theoretical predictions (curves) based on the linearised analytical scaling law of (4.1). The strong agreement of simulation data with the scaling law implies that, for all parameters employed here, the maximum velocity is already reached in the quasi-linear phase of the blob evolution.

Figure 2. Maximum radial blob velocity $v_{\max }^x$ obtained from nonlinear full-F multi-species simulations as a function of the relative density ratio $a$ and mass ratio $\mu$, in comparison with the approximate analytical scaling law.

This scenario of blob propagation is illustrated in figure 3. For $\tau _i = \tau _j = 0$ the evolution is up–down mirror symmetric with respect to the ‘poloidal’ coordinate $y$ in the $(x,y)$ plane perpendicular to $\boldsymbol {B}$. This redundancy is here used to illustrate the differences in evolution for two pairs of impurity concentration and mass. The $(x,y)$ plane with $L_x = L_y = 500 \rho _0$ is divided into an upper half ($y > 250$) that shows the case for the $\mu _j=2$, $a_j=0.001$ parameters, and a lower half ($y<250$) that shows the case for the $\mu _j=20$, $a_j=0.5$ parameter set. A snapshot of the electron density $n_e(x,y)$ is displayed for three different simulation times ($t=0$, $t=983$ and $t=1224$ in units of $\varOmega _i^{-1}$). The blob propagation for the case with high mass-to-charge density $\mu _j$ is visibly slower, but both blobs have for the later times retained the characteristic laminar mushroom top-like shapes. The overall scaling of the full-F blob velocity is a generalisation of the isotopic mass scaling found in delta-F simulations in Meyer & Kendl (Reference Meyer and Kendl2017) towards non-isotopic arbitrary impurities.

Figure 3. Propagation of two cold ($\tau _i = \tau _j = 0$) full-F blobs containing impurities with parameters ($\mu _j=2, a_j=0.001$) and ($\mu _j=20, a_j=0.5$) at the initial time (left) and in two evolved phases. The upper halves (of otherwise up–down symmetric blobs) show the density $n_e(x,y)$ for the first parameter case, and the lower halves for the second case with higher effective mass density, which clearly slows the propagation. The linear density colour scale ranges from $N_{e0} = 1$ (light/grey) to $(N_{e0}+\Delta N)= 2$ (dark/blue). A ($150\times 150$) $\rho _0$ central section of the ($500\times 500$) $\rho _0$ simulation domain is here shown.

The present isothermal full-F gyro-fluid model also allows simulations with finite ion and impurity temperatures given by the parameters $\tau _i$ and $\tau _j$, and the associated FLR effects. Then the poloidal up–down symmetry is broken, and quantitative comparisons with maximum velocity scaling laws are less exact. A complete depiction of FLR effects is also only feasible with the full-F full-k gyro-fluid model, for which a multi-species code implementation is not yet available. Further, a final goal of seeded single blob studies is the understanding of turbulent filamentary transport in the SOL of fusion plasmas, which can be seen as an ensemble of blobs of a variety of sizes that are generated dynamically around the separatrix regions. For these reasons, here only an exemplary demonstration of a warm full-F impurity enriched plasma blob is presented. In later studies, impurity transport will be studies in edge–SOL coupled 3-D full-F full-k turbulence codes, including both the drift wave driven closed-field-line edge region and the open-field-line (‘blobby’) SOL region.

In figure 4 the evolution of a blob perturbation with impurity concentration $a_j = 0.01$, relative mass/charge ratio $\mu _j = 2$ and relative temperatures $\tau _i = \tau _j = 2$ is shown at three different stages of the dynamical evolution. In addition to the density $n_e(x,y)$, contour lines of the asymmetric dipolar electric potential $\phi (x,t)$ are also drawn in the figure. The overall structure and dynamics are similar to those of warm blobs composed of only a single ion species.

Figure 4. Propagation of a warm full-F blob containing impurities with concentration $a_j=0.01$, $\mu _j=2$ and finite temperature ratios $\tau _i = \tau _j = 2$. In addition to the colour scale plot of density $n_e(x,y)$, also exemplary isocontour lines of the electric potential $\phi$ are shown.

5. Blob interaction with a non-trace impurity cloud

The former simulations presented above have assumed an equally distributed background concentration of impurity ions. The ability of the multi-species full-F gyro-fluid model to also treat self-consistent interactions of the main plasma dynamics with localised non-trace impurity is now demonstrated in simulation case ‘A’: on a background plasma with homogeneous $N_e = N_i = 1$ and $\tau _i = 1$, an impurity cloud is added which is localised to a radial region around $x_0 = 250 = 0.5 L_x$ with a Gaussian width of $\sigma _{{\rm imp}} = 10$, and extended in the $y$ direction. The impurity background concentration is set to $a_j = 0.001$, $\mu _j=2$ and $\tau _j =1$. The maximum amplitude of the initial impurity distribution is $\Delta N_j / N_{j0} = 5$.

In addition, a warm ($\tau _i=1$) quasi-neutral Gaussian blob composed only of main ion species density and electrons is added at $(x,y) = (0.4 L_x, 0.5 L_y)$ with width $\sigma = 10$ and relative amplitude $\Delta N_i/N_{i0}=1$, which again propagates radially outwards (to larger $x$), and into the impurity cloud. The self-consistent evolution is shown in figure 5, where the impurity cloud density $N_j(x,y)$ is depicted in (red) shades, whereas the main plasma blob density $N_i(x,y)$ is visualised by (blue) isocontour lines.

Figure 5. Case ‘A’ (for parameters see main text): illustration of a self-consistent interaction of a main plasma species warm blob (blue isocontours) with a localised low concentration cloud of non-trace impurities (red). Temperature ratios are $\tau _i = \tau _j = 1$.

The advection by the combined E×B velocity through the electric potential computed by the multi-species polarisation equation leads to transport of the impurity density around the main plasma blob. The initial symmetry of the impurity cloud $N_{j0}(x)$ in the $y$ direction is broken by the impact, and a local deformation occurs. Mixing during this quasi-linear laminar phase of blob evolution occurs here only in the trailing dipolar vortex lobes. The overall impurity cloud is shifted (‘pinched’) inwards by the passing blob. This evolution for exemplary low concentration $a_j = 0.001$ and mass $\mu _j=2$ is still highly similar to passive trace impurity advection, and the trajectory of the main blob is not noticeably influenced.

A perhaps at first glance surprising or counter-intuitive effect can be observed with higher relative density of the localised impurity wall on an impinging main species blob. The results shown in figure 2 imply a slower maximum blob velocity in the presence of impurities. However, in that former case of a background distributed impurity the higher inertia entered through an increased effective mass of the mixed (impurity and main ion) blob perturbation. In contrast, in the present scenario, the initial blob perturbation is composed purely of main ions and enters the initially poloidally symmetric cloud of impurities at some time during its propagation.

As can be seen from figure 5, the impurity density is advected along iso-contours of the blob electric potential (as a streamfunction for the E×B velocity), which does not influence the effective mass of the main blob. But once the impurity wall is perturbed by the impinging blob, the developing poloidal impurity asymmetry gives rise to additional interchange driven baroclinic vorticity generation. The primary blob potential and the secondary impurity asymmetry potential superimpose, and together act on both blob and impurity advection.

The effect on the main blob propagation is demonstrated by a series of simulations with varying relative impurity wall concentration and location. All simulations are initialised with a cold ($\tau _i=0$) main ion blob of Gaussian width $\sigma _i = 10$ in units of $\rho _s$ and amplitude $A_i=1$ relative to the constant main ion background density.

The simulation case ‘B’ now initialises the cold ($\tau _j=0$) impurity cloud on the radial position $x_{0j} = 0.45 L_x$ outwards of the main blob ($x_{0i} = 0.3 L_x$) of same Gaussian width $\sigma _j = 10$, and with a homogeneous initial impurity distribution in the (poloidal) $y$ direction. The amplitude of the cloud is $A_j = 5$ relative to a constant impurity background with concentration $a_j=0.1$, so two orders of magnitude higher than for the quasi-trace case ‘A’ and similar in magnitude to the main ion density. The larger initial separation of $\Delta x = 0.15 L_x$ compared with case ‘A’ (with $0.1 L_x$) ensures that the blob accelerates mostly before it reaches the cloud and reaches its maximum velocity around the left, increasing slope of the cloud. The maximum centre of mass velocity blob in this impurity cloud impact case ‘B’ is $v_{\max }^B = 0.0416$, compared with the velocity $v_{\max }^0 = 0.0532$ of the pure blob without any impurity cloud. This reduced velocity is comparable to those expected in the background mixture scenario depicted in figure 2 for values of $\mu _j=3$ and a background concentration in the (here in case ‘B’ inhomogeneous) range $0.1 \leqslant a_j \leqslant 0.5$.

The interpretation of the blob slowing down is here, however, not based on the effective mass of the initial blob (which is purely main species ions), but on the counteracting inverted dipole potential added at the main ion blob front by the impurity ‘hole’, which is indented into the impurity cloud, as illustrated by the left panel ($t=710$) in figure 5. The total accelerating dipole field is therefore reduced compared with the pure blob without impurities. Again, a significant effect is only expected for rather high impurity concentrations, similar to the mixed background impurity cases above.

For the next series of simulations (‘case C’, see figure 6), the cold ($\tau _j=0$) impurity cloud is initialised on the same radial position $x_{0j} = 0.3 L_x$ as the main blob and with same Gaussian width $\sigma _j = 10$. This implies that the blob is accelerated to the right over the decreasing impurity density slope. The amplitude of the cloud is $A_j = 5$ relative to a constant impurity background with concentration $a_j$ varied between $0.03$ and $0.2$, relative to the main ion density. The impurity mass $\mu _j$ is varied between 1 and 7. For these initialisations, the accelerating blob mostly pushes impurity density out to the right (in its propagation direction), which again leads to a poloidal impurity density asymmetry, but now in a blob-like form near the front of the main ion blob. Both dipole fields are superimposed constructively and add to the acceleration of the primary main ion blob. This situation is graphically illustrated in figure 7.

Figure 6. Case ‘C’: maximum velocity $v_{\max }^x$ of a blob started inside a radially localised impurity cloud. The blob initially propagates down the impurity gradient, which for light impurities (low $\mu _j$) adds constructively to the interchange acceleration of the blob. (Lines are interpolation between simulation data points to guide the eye.)

Figure 7. (a,b) In case ‘B’ (and similar to case ‘A’) the main ion blob propagates towards the impurity wall. The $E\times B$ advection along isocontours of the dipolar potential (here shown by blue and red contours) pushes impurities out to the right. This is similar to as if an impurity blob appeared on the right side of the wall and a hole on the left side, each of which generates its own electric field with opposite signs of the dipolar potentials. The impurity blob-like perturbation always is ahead of the main ion blob, which coincides more closely with the impurity hole: both main blob and impurity hole potentials are destructively superimposed, so that the main blob speed is reduced. (c,d) In case ‘C’ the main ion blob is initialised inside the impurity wall. Now the arising impurity blob on the right side of the wall mostly coincides with the main blob, whose potentials constructively superimpose and add to the propagation velocity.

This leads to the results depicted in figure 6: for light (small effective mass $\mu _j$) and low to medium impurity density, the main ion blob achieves a higher maximum velocity compared with the case without any impurities (corresponding to the points $a_j=0$ in the plot), and outward transport is locally enhanced. Again, the effect is rather negligible for low concentrations ($a_j \leqslant 0.01$), but could be effective, for example, for gas puffing in the SOL of fusion experiments. It should here be remarked that, in case of the presence of neutral components, the effects of ionisation and collisions with neutrals would also affect the blob propagation. These are not included in the present model.

An acceleration of the main species blobs by an impurity cloud may appear counter-intuitive at first glance, but is here caused by the locally co-aligned impurity and main ion density gradients. Such situations can appear in an experiment when the maximum impurity density is achieved around the separatrix region of a fusion experiment (by inward pinching), where, again, the major blob production is also expected to occur.

Later stages of blob development can lead to turbulent regions, for example generated by Kelvin–Helmholtz unstable blob fronts, or by additional drift wave effects when parallel dynamics is included in the model. For such small-scale dynamics the full treatment of FLR effects is desirable for consistency, therefore such studies have to be postponed until future availability of full-F full-k codes.

The present simulations show qualitative similarity to the results of the full particle-based PIC computations from Hasegawa & Ishiguro (Reference Hasegawa and Ishiguro2017), which used different parameters (such as an artificial low ion to electron mass ratio of 100/1 in the PIC model, and small Larmor radii). A similar advection dynamics is, however, also very visible in the (by nature noisy) PIC simulations, which might in future (with extended PIC simulations) possibly allow for cross-verification between full particle and gyro-fluid or gyro-kinetic models.