Global turbulence simulations of the tokamak edge region with GRILLIX

TL;DR

This paper extends the GRILLIX fluid-turbulence code to a global model for the tokamak edge, capturing large fluctuations and complex geometries, and compares turbulence characteristics between circular and diverted geometries.

Contribution

The authors develop a global, electromagnetic turbulence model based on the drift-reduced Braginskii equations without Boussinesq approximation, incorporating realistic geometries using FCI techniques.

Findings

Turbulence is highly intermittent with fluctuations up to 40%.

Diverted geometry reduces fluctuation levels compared to circular geometry.

The Boussinesq approximation has minimal impact on turbulence results.

Abstract

Turbulent dynamics in the scrape-off layer (SOL) of magnetic fusion devices is intermittent with large fluctuations in density and pressure. Therefore, a model is required that allows perturbations of similar or even larger magnitude to the time-averaged background value. The fluid-turbulence code GRILLIX is extended to such a global model, which consistently accounts for large variation in plasma parameters. Derived from the drift reduced Braginskii equations, the new GRILLIX model includes electromagnetic and electron-thermal dynamics, retains global parametric dependencies and the Boussinesq approximation is not applied. The penalisation technique is combined with the flux-coordinate independent (FCI) approach [F. Hariri and M. Ottaviani, Comput.Phys.Commun. 184:2419, (2013); A. Stegmeir et al., Comput.Phys.Commun. 198:139, (2016)], which allows to study realistic diverted geometries with X-point(s) and general boundary contours. We characterise results from turbulence simulations and investigate the effect of geometry by comparing simulations in circular geometry with toroidal limiter against realistic diverted geometry at otherwise comparable parameters. Turbulence is found to be intermittent with relative fluctuation levels of up to 40% showing that a global description is indeed important. At the same time via direct comparison, we find that the Boussinesq approximation has only a small quantitative impact in a turbulent environment. In comparison to circular geometry the fluctuations are reduced in diverted geometry, which is related to a different zonal flow structure. Moreover, the fluctuation level has a more complex spatial distribution in diverted geometry. Due to local magnetic shear, which differs fundamentally in circular and diverted geometry, turbulent structures become strongly distorted in the perpendicular direction and are eventually damped away towards the X-point.

Tables

Table 1. Table 1: Parameters for analytic MMS functions used as inputs for different dynamical fields.

	$k_{\rho }$	$(k_{\theta },\delta \theta )$	$(k_{\phi },\delta \phi )$	$\omega$
${\theta }_n$	$1$	$(1,\mathrm{\hspace{0.17em}0})$	$(1,\mathrm{\hspace{0.17em}0})$	$100$
${\xi }_e$	$1$	$(1,\mathrm{\hspace{0.17em}1.5})$	$(1,\mathrm{\hspace{0.17em}0.5})$	$73$
$\varphi$	$2$	$(2,\mathrm{\hspace{0.17em}0})$	$(1,\mathrm{\hspace{0.17em}0})$	$80$
$u_{\parallel }$	$3$	$(1,\mathrm{\hspace{0.17em}0})$	$(1,\mathrm{\hspace{0.17em}0})$	$65$
$A_{\parallel }$	$2$	$(1,\mathrm{\hspace{0.17em}0})$	$(1,\mathrm{\hspace{0.17em}0})$	$88$

Table 2. Table 2: Numerical parameters used for convergence check. w θ subscript 𝑤 𝜃 w_{\theta} is poloidal decay length in radians and w ρ subscript 𝑤 𝜌 w_{\rho} radial decay length in units of R 0 subscript 𝑅 0 R_{0} for penalisation function, which is parametrized via tanh \tanh functions.

	$h_{⟂}$	$N_{pol}$	$\nu$	$\mu$	$w_{\theta }$	$w_{\rho }$
coarse	$2.0$	$16$	$500.0$	$5.0\cdot {10}^{-2}$	$0.1$	$6.0\cdot {10}^{-3}$
nominal	$1.0$	$32$	$10.0$	$2.5\cdot {10}^{-2}$	$7.5\cdot {10}^{-2}$	$4.0\cdot {10}^{-3}$
fine	$0.67$	$48$	$2.0$	$1.5\cdot {10}^{-2}$	$5.0\cdot {10}^{-2}$	$3.0\cdot {10}^{-3}$

Equations

$$\frac{d_u}{dt}{\theta }_n+\nabla \cdot \left(\mathbf{b}{u}_{\parallel }\right)=$$

$$+\frac{1}{n}\left[{\mathcal{D}}_n\left(n\right)+S_n\right],$$

$$\nabla \cdot \left[\frac{n}{B^2}\frac{d_u}{dt}{\nabla }_{\perp }\varphi \right]=$$

$$+{\mathcal{D}}_w\left(\Omega \right),$$

$$\frac{d_u}{dt}{u}_{\parallel }=$$

$$\frac{\partial }{\partial t}{\psi }_{\parallel }+\mu \left({\mathbf{v}}_E\cdot \nabla +v_{\parallel }{\nabla }_{\parallel }\right)\left(\frac{j_{\parallel }}{n}\right)=$$

$$+T_e\left[{\nabla }_{\parallel }{\theta }_n+1.71{\nabla }_{\parallel }{\xi }_e\right]-{\nabla }_{\parallel }\varphi ,$$

$$+{\mathcal{D}}_p\left({\psi }_{\parallel }\right)$$

$$\frac{d_v}{dt}{\xi }_e=$$

$$+\frac{2}{3}0.71\frac{1}{n}\nabla \cdot \left(\mathbf{b}{j}_{\parallel }\right)$$

$$\frac{1}{T_e}\left[{\mathcal{D}}_t\left(T_e\right)+S_t\right],$$

$${\nabla }_{\perp }^2{A}_{\parallel }=-j_{\parallel }.$$

$${\mathbf{u}}_E\cdot \nabla f=$$

$$\mathcal{C}(f)=$$

$${\mathcal{D}}_f:={\nu }_{\perp ,f}{\nabla }_{\perp }^{2N}f+{\nu }_{\parallel ,f}\nabla \cdot \left(\mathbf{b}{\nabla }_{\parallel }f\right),$$

$$u_{\parallel }\gtrless \pm \sqrt{T_e},$$

$$j_{\parallel }=0,$$

$$\varphi =\Lambda T_e$$

$$-{\chi }_{\parallel }{\nabla }_{\parallel }{T}_e={\gamma }_e{T}_e{u}_{\parallel },$$

$${\nabla }_{\parallel }n=0,$$

$${\nabla }_{\parallel }{f}_{i,k+\frac{1}{2},j}:=\frac{f_{k+1}\left({\mathbf{\gamma }}_{i,j}\left(\frac{\Delta \phi }{2}\right)\right)-f_k\left({\mathbf{\gamma }}_{i,j}\left(-\frac{\Delta \phi }{2}\right)\right)}{s_{i,j}\left(\frac{\Delta \phi }{2}\right)+s_{i,j}\left(-\frac{\Delta \phi }{2}\right)},$$

$$\frac{d{\mathbf{\gamma }}_{i,j}}{d\phi }=$$

$$\frac{d{s}_{i,j}}{d\phi }=$$

$$\frac{\partial }{\partial t}f=\left(1-\chi \right)F_f+\frac{\chi }{ϵ}\left(f_P-f\right),$$

$${{\nabla }_{\parallel }f|}_{sheath}=\alpha .$$

$$f_P=\{\begin{array}{ll}\left|\zeta \right|\left(f^{-}+s^{-}\alpha \right)+(1-\left|\zeta \right|)\frac{f^{+}+f^{-}}{2}& \text{for: }\zeta \ge 0,\\ \left|\zeta \right|\left(f^{+}-s^{+}\alpha \right)+(1-\left|\zeta \right|)\frac{f^{+}+f^{-}}{2}& \text{for: }\zeta <0.\end{array}$$

$$\frac{\chi }{ϵ}{\varphi }^{t+1}-\left(1-\chi \right)\nabla \cdot \left(\frac{n^{t+1}}{B^2}{\nabla }_{\perp }{\varphi }^{t+1}\right)=\frac{\chi }{ϵ}\left(\Lambda T_e^{t+1}\right)-\left(1-\chi \right){\Omega }^{t+1},$$

$${\beta }_0{A}_{\parallel }^{t+1}-\frac{\mu }{n^{t+1}}{\nabla }_{\perp }^2{A}_{\parallel }^{t+1}={\psi }_{\parallel }^{t+1},$$

$$f^{t+1}\left(11+6\Delta t\frac{\chi }{ϵ}\right)=\sum _{i=0\dots 2}{a}_i{f}^{t-i}+b_i\Delta t\left[\left(1-\frac{\chi }{ϵ}\right)F^{t-i}(f^{t-i})+\frac{\chi }{ϵ}{f}_P^{t-i}\right],$$

$$c_{1}f-c_2\nabla \cdot \left(c_3{\nabla }_{\perp }f\right)=b,$$

$$\frac{{\partial }^2}{\partial t^2}\left(\tilde{n}-\tilde{\varphi }\right)=v_{SAW}^2{\nabla }_{\parallel }^2\left(\tilde{n}-\tilde{\varphi }\right),$$

$$\frac{\partial }{\partial t}\varphi =\frac{1}{k_{\perp }^2{\eta }_{\parallel }}\nabla \cdot \left(\mathbf{b}{\nabla }_{\parallel }\varphi \right).$$

$$R_0=$$

$$n_0=$$

$${\beta }_0=$$

$${\chi }_{\parallel 0}=$$

$$S_n=c_n\exp \left(-\frac{(\rho -{\rho }_{src}{)}^2}{w_{src}^2}\right)\left({⟨n⟩}_{\rho }-n_{target}\right)$$

Full text

††journal:  

1 Introduction

Understanding the complex multi-physics of the edge region – the scrape-off layer (SOL) and closed-field line region immediately near the separatrix – is of critical importance for the development of fusion energy. Due to the relative stiffness of the internal profiles of temperature and density, the core values and therefore the overall fusion performance is strongly determined by the edge profiles. Furthermore, in a fusion reactor there will be a large exhaust of particles and heat due to imperfect confinement. This plasma exhaust is directed towards divertor target plates, and preventing these fluxes from exceeding engineering limits, above which the performance and lifetime of the reactor is significantly reduced, is a high-priority area of fusion research. Prediction of these heat fluxes for future devices such as ITER or DEMO is complicated by uncertainty in extrapolation of the width of the exhaust channel from current devices [1, 2].

Modelling of the edge plasma is a significant challenge due to the highly-coupled interplay of multiple different physics regimes and disparate spatial and temporal scales. Magnetized plasma physics, complex magnetic geometry, neutral physics and momentum transfer, atomic and molecular chemistry, radiation from excited states, wall recombination, surface chemistry and impurity sputtering all can affect the edge plasma. Furthermore, the edge can exhibit phenomena over a large range of spatial and temporal scales – from the formation of small-scale intermittent turbulent filaments to large-scale long-timescale effects such as equilibration of the background in response to the magnetic and wall geometry. Inclusion of an extended physics set or finer spatial and temporal scales typically improves the accuracy of the code with respect to experiment, but at the expense of increased computational cost. Within the subset of codes based on the multi-fluid approximation, the two broad classes are ‘transport’ and ‘turbulence’ codes. Transport codes such as SOLPS code [3] include a significant range of multi-physics but do not treat turbulent transport self-consistently. Instead, the effects of turbulence are approximated via an effective diffusion, which remains an ad-hoc input. In contrast, turbulence codes self-consistently treat turbulence by evolving the 3D Braginskii models [4] – at the cost of increased runtime and/or a reduced physics set.

Several recent projects aim at developing fluid-turbulence codes, of which we note the GBS [5, 6, 7], HERMES (BOUT++) [8], TOKAM3X [9], GDB [10] and GRILLIX [11] projects. In contrast to the other codes mentioned, the GRILLIX project is notable for its use of the flux-coordinate independent (FCI) approach [12, 13, 14]. This method prevents the issue of coordinate singularities at the separatrix and X-point which arise from the use of field- or flux-aligned coordinates. GRILLIX employs a cylindric grid $(R_{i},\varphi_{k},Z_{j})$ where parallel operators are discretised via field line tracing between toroidal planes and field line map interpolation within each plane. This allows for the use of a single consistent method to be used for all grid points, including the possibility of (possibly multiple) X-points. Furthermore, the use of Cartesian grids prevents resolution imbalances between the outboard mid plane and the X-point region, allowing the dynamics around the X-point to be investigated with high accuracy. To allow for the treatment of general non-conformal boundaries, the penalization method is used to enforce the desired boundary conditions.

In this paper, the extension of GRILLIX by electromagnetic and electron thermal dynamics is presented, with the resulting model being a global drift reduced Braginskii model. Here ’global’ means that parametric dependencies are kept and that the Boussinesq approximation is not applied, i.e. nowhere a separation is made between fluctuations and background for the density and temperature. This is needed for the consistent description of high amplitude fluctuations which are regularly observed in experiments [15]. Recently several codes were adapted to abolish the Boussinesq approximation also for turbulence applications [6, 10], whereat now the geometrical complexity of diverted equilibria is additionally introduced in GRILLIX. The Boussinesq approximation has often been studied at isolated sub problems, i.e. blob propagation [16, 17, 18, 19]. Here we investigate its impact in a fully turbulent environment and find that it has only minor quantitative impact on results. Using parameters characteristic for the COMPASS tokamak [20] turbulence simulations are carried out with GRILLIX. The turbulence features intermittency and exhibits large relative fluctuation levels, which shows that a global description is generally important. The impact of geometry is studied by comparing simulations in circular geometry with toroidal limiter against realistic diverted geometry at otherwise comparable parameters. The fluctuation level is reduced in diverted geometry owed to a different zonal flow structure. Moreover, the fluctuation level exhibits a more complex spatial distribution, which we explain as a consequence of local magnetic shear, which fundamentally differs in the edge region between circular and diverted geometry. Strong local magnetic shear causes a distortion of turbulent structures, which become subsequently subject to strong perpendicular dissipation. The X-point, where magnetic shear becomes locally very strong, thereby tends to disconnect the low field side from the high field side, where curvature acts as stabilizing [21]. Therefore, stronger poloidal asymmetries in the fluctuation level are observed in diverted geometry.

The remainder of this paper is organized as follows: In section 2 we present the physical model employed in GRILLIX, a global 3D drift reduced Braginskii model [22]. With respect to the previous version [11] GRILLIX has been extended by electromagnetic dynamics, electron temperature dynamics and the Boussinesq approximation has been relaxed, which enables to simulate plasma turbulence globally, i.e. without splitting quantities into background and fluctuations. The implementation of important new features is described in section 3. A geometric multigrid solver for the generalised 2D perpendicular Helmholtz equation allows an efficient treatment of the new electromagnetic terms and relaxation of the Boussinesq approximation. The treatment of sheath boundary conditions at the divertor/limiter plates via penalization techniques was motivated from the GDB code [10, 23] and we give a generalisation to diverted geometries. The extended model and new features are verified by analytic means and the method of manufactured solutions (MMS) [24]. In section 4 we characterize edge turbulence with GRILLIX simulations and clarify the impact of geometry by comparing simulations in circular geometry with toroidal limiter against simulations in diverted geometry at otherwise comparable parameters. A summary and outlook is given in section 5.

2 Physical model

2.1 Global drift reduced Braginskii equations

Based on the assumptions of short mean free paths, i.e. $\lambda_{c}\ll R_{0}$, the drift reduced Braginskii model describes plasma dynamics of low frequency ($\omega\ll\Omega_{i}$) in comparison to the ion cyclotron frequency $\Omega_{i}$ [4, 22].It is suitable to describe turbulence at low temperature in the edge region self consistently. As further practical approximations cold ions ($T_{i}\ll T_{e}$) are assumed in GRILLIX, and whereas electromagnetic effects are kept in Ohm’s law, magnetic flutter is neglected, i.e. transport is assumed to be electrostatic and the description of certain phenomena like Edge Localised Modes (ELM’s) is excluded.

The following normalisation is employed: Time $t$ is normalised against $R_{0}/c_{s0}$ with the sound speed $c_{s0}:=\sqrt{T_{0}/M_{i}}$ at some reference temperature $T_{0}$. Parallel scales $x_{\parallel}$ are normalised against major radius $R_{0}$ and perpendicular scales $x_{\perp}$ against the sound Larmor radius $\rho_{s0}:=c\sqrt{T_{0}M_{i}}/(eB_{0})$ with $B_{0}$ the magnetic field strength on axis and $M_{i}$ is the ion mass. The dynamical variables in GRILLIX are density $n$ normalized against some reference density $n_{0}$, electron temperature $T_{e}$ against $T_{0}$, parallel ion $u_{\parallel}$ and electron $v_{\parallel}$ velocities against $c_{s0}$, parallel current $j_{\parallel}$ against $en_{0}c_{s0}$, electrostatic potential $\phi$ against $T_{0}/e$ and the parallel component of the electromagnetic potential $A_{\parallel}$ against $\beta_{0}B_{0}\rho_{s0}$ with $\beta_{0}:=4\pi n_{0}T_{0}/B_{0}^{2}$ the dynamical plasma beta at reference values. In order to preserve its positivity the logarithms of normalised density $\theta_{n}:=\log n$ and temperature $\xi_{e}:=\log T_{e}$ are evolved in time. Finally, the normalised set of equations implemented in GRILLIX are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Equations (1-6) are the electron continuity equation, vorticity equation or quasineutrality condition, parallel momentum equation, Ohm’s law, electron temperature equation and Faraday’s law respectively. The total time derivatives contain advection by the E$\times$B drift and the parallel velocity, i.e. $\frac{d_{u}}{dt}:=\frac{\partial}{\partial t}+\mathbf{u}_{E}\cdot\nabla+u_{\parallel}\nabla_{\parallel}$ and $\frac{d_{v}}{dt}:=\frac{\partial}{\partial t}+\mathbf{u}_{E}\cdot\nabla+v_{\parallel}\nabla_{\parallel}$. The parallel gradient $\nabla_{\parallel}=\mathbf{b}\cdot\nabla$ is taken along the static equilibrium magnetic field $\mathbf{B}$ with unit vector $\mathbf{b}:=\mathbf{B}/B$, i.e. magnetic flutter is not included. As auxiliary variables the generalised vorticity $\Omega:=\nabla\cdot\left(\frac{n}{B^{2}}\nabla_{\perp}\phi\right)$ and the generalised parallel electromagnetic potential $\psi_{\parallel}:=\beta_{0}A_{\parallel}+\mu\frac{j_{\parallel}}{n}$ have been introduced. Within the FCI approach (see section 3.1) the advection by E$\times$B velocity and the curvature operator can be written as:

[TABLE]

with the Jacobi bracket $\left[\phi,f\right]_{R,Z}:=\partial_{R}\phi\partial_{Z}f-\partial_{Z}\phi\partial_{R}f$. The dimensionless parameters of the system are the drift scale $\delta:=R_{0}/\rho_{s0}$, electron to ion mass ratio $\mu:=\frac{m_{e}}{M_{i}}$, the dynamical plasma beta $\beta_{0}$, reference parallel resistivity $\eta_{\parallel 0}:=0.51\mu/(\tau_{e0}(c_{s0}/R_{0}))$ with $\tau_{e0}$ the electron-ion collision time at reference values and the reference parallel heat conductivity $\chi_{\parallel 0}:=3.15/\mu\tau_{e0}(c_{s0}/R_{0})$. Particle and thermal source terms $S_{n},S_{t}$ have been added, and mostly for numerical reasons also dissipation terms:

[TABLE]

with constant coefficients $\nu_{\perp,f},\,\mu_{\parallel,f}$, and $N$ controlling the order of perpendicular dissipation. Note that in order to ensure conservation of particles the dissipation in the continuity equation (1) does not act on the logarithm of the density $\theta_{n}$, but the density $n$.

The model is global in the sense that no separation of variables between a background part and fluctuating part is made and the dependency of the parallel resistivity and heat conduction on density and temperature is kept. Moreover the dependency on the density in the polarization term of the vorticity equation (2) is also kept, i.e. the Boussinesq approximation is not applied. The model conserves energy apart from the fact that we have neglected advection by the polarization velocity, which has been shown to have only a minor effect on conservation of energy [25].

2.2 Boundary conditions

The simulation domain in GRILLIX is usually bounded by an inner (core) limiting flux surface, an outer (wall) limiting flux surface and depending on geometry by limiter or divertor baffle plates.

Sheath physics determining the boundary conditions for the divertor/limiter is a rich topic by itself. A sophisticated set of boundary conditions is given e.g. in [26], where an inclination of the magnetic field with respect to the target plates is taken into account, but more commonly just Bohm boundary conditions [27] are employed. Moreover, the treatment of sheath boundaries within the FCI approach is also numerically cumbersome for which a penalization method is employed in GRILLIX described in section 3.2. Due to all these complexities we restrict ourselves firstly to the relatively simple but robust set of insulating sheath boundary conditions. This assumption implies, e.g. that blobs propagate purely according to the inertial scaling $v_{blob}\propto w_{blob}^{1/2}$, with $v_{blob}$ the radial blob velocity and $w_{blob}$ the blob width, whereas with Bohm boundary conditions larger blobs would propagate slower according to the sheath connected scaling $v_{blob}\propto w_{blob}^{-2}$ [28]. Therefore, the assumption of an insulating sheath might appear strong, but is of relevance for modelling detached conditions where resistivity in the front of the target plates is strongly enhanced due to the presence of neutrals.

Finally, the boundary conditions employed in GRILLIX are:

[TABLE]

where the upper/lower sign denotes if the direction of magnetic field is directed towards/away from target plates, $\Lambda\approx 0.5\ln\left(\frac{M_{i}}{2\pi m_{e}}\right)$ is the sheath floating potential and $\gamma_{e}\approx 2.5$ the effective electron sheath transmission factor. As the continuity equation is of hyperbolic nature in the parallel direction the boundary condition on the density should be as unrestricted as possible, and we found that a homogeneous Neumann boundary condition is numerically more robust than extrapolation $(\nabla_{\parallel}^{2}n=0)$. The generalised vorticity and electromagnetic potential are obtained consistently with the electrostatic potential respectively the parallel current as described in section 3.2.

In the radial direction either homogeneous Neumann boundary conditions (for $n$, $T_{e}$, $u_{\parallel}$) are applied or homogeneous Dirichlet boundary conditions (for $A_{\parallel}$, $j_{\parallel}$, $\Omega$). An exception is the electrostatic potential which is set at the wall to $\left.\phi\right|_{wall}=\Lambda T_{e}$. In order to avoid fluxes of energy and particles due to $E\times B$ drifts through the core, the potential has to be constant on the inner limiting flux surface, and we set it to $\phi_{core}=\Lambda\left\langle T_{e}\right\rangle_{LCFS}$, where $\left\langle T_{e}\right\rangle_{LCFS}$ is the zonal averaged electron temperature on the last closed flux surface. The motivation for this stems from the fact that due to the sheath boundary conditions the potential follows roughly $\phi\sim\Lambda T_{e}$ in the SOL, and the chosen boundary condition does therefore not allow a global radial electric field in the closed flux surface region. The large scale radial electric field in the closed field line region is determined by effects that are not yet included in the GRILLIX model, e.g. ion pressure gradient [29, 30]. Therefore, the self-consistent modelling of the global radial electric field in the closed field line region is postponed until at least ion temperature effects will be taken into account in GRILLIX. We want to note that from the geometrical point of view GRILLIX is able to deal easily with the full tokamak including the core region with O-point [11]. Whereas this would exclude a possibly spurious influence from core boundaries and therefore allow a more self-consistent approach, the many additional grid points would pose a large computational overhead.

3 GRILLIX

3.1 Spatial discretisation

GRILLIX is based on the FCI approach [12, 13, 14] in a toroidally staggered framework which is described in detail in [11] and is therefore here only reviewed very shortly. For tokamak geometries a cylindric grid $(R_{i},\varphi_{k},Z_{j})$ is employed being Cartesian within poloidal planes. Based on the assumption of strong toroidal field ($B_{pol}/B_{tor}\ll 1$), the discretisation of perpendicular operators remains within poloidal planes for which second order finite difference methods are used. The Jacobi bracket is discretised according to the Arakawa scheme [31] and the discretisation of the non-linear polarization term in eq. (2) is described in [25]. The discrete parallel gradient is computed at toroidally staggered positions $\varphi_{k+\frac{1}{2}}$ according to a finite difference along magnetic field lines (see fig. 1):

[TABLE]

where $\boldsymbol{\gamma}_{i,j}(\varphi)$ is the poloidal projection of the characteristic along field line and $s_{i,j}(\varphi)$ the associated length along field line defined as the solution of the following ordinary differential equations which are solved via a Runge-Kutta integrator [32]:

[TABLE]

i.e. corresponding map points are computed by tracing along magnetic field lines. The values on map points are obtained from a 3rd order bi-polynomial interpolation within the poloidal planes $\varphi_{k}$ and $\varphi_{k+1}$. We note that the magnetic field is assumed axisymmetric in GRILLIX, which is however not a general constraint for the FCI approach [33, 34]. In the same spirit operators are established that map quantities between the grid and the staggered dual grid and vice versa. The parallel divergence operator $\nabla\cdot\left(\mathbf{b}f\right)$ is obtained in its discrete version via the support operator method [35, 36] as described in [11]. The structure of the equations suggest that $n,\theta_{n},\xi_{e},T_{e},\phi$ and $\Omega$ are co-located to the canonical grid whereas $u_{\parallel},v_{\parallel},j_{\parallel},A_{\parallel}$ and $\psi_{\parallel}$ are co-located to the staggered grid.

3.2 Penalization for sheath boundary conditions

The boundaries at the sheath are in general neither conformal with the grid nor aligned with the exceptional parallel direction along the magnetic field line, which makes their treatment numerically difficult and cumbersome. In such situations penalization techniques have proven themselves also for plasma fluid codes [37, 38]. A combination of the FCI with penalization was firstly employed in GDB for limited plasmas [10, 23] and we implemented in GRILLIX a generalization allowing to deal also with diverted plasmas.

Equations (1-5) are each modified according to:

[TABLE]

where $f$ represents here the dynamical variables respectively and $F_{f}$ the corresponding terms according to the Braginskii model. $\chi$ is a characteristic function, which is [math] in the physical domain and $1$ in the boundary region, where we choose in practice a smooth transition [37] across the boundary based on $\tanh$ functions (see fig. 5 left column for examples). $\epsilon\ll 1$ is the penalization parameter such that in the region where $\chi\approx 0$ eq. (13) approximates the original physical equation, whereas in the region $\chi\approx 1$ the variable $f$ is strongly damped to a prescribed function $f_{P}$. Via suitable choice for $f_{P}$ different boundary conditions can be realized.

As an illustrative example we discuss here our implementation for general Neumann boundary conditions, i.e.:

[TABLE]

Firstly, we define an additional function $\zeta$, which is $1$ in the penalization region where the magnetic field is pointed towards the target and $-1$, where the magnetic field is pointed away from the target with possibly, i.e. in toroidal limiter geometry, a smooth transition between both regions (see fig. 5 right column for examples). Secondly, we denote for some grid point $f^{\pm}:=f_{k\pm 1}\left(\boldsymbol{\gamma}_{i,j}(\pm\Delta\varphi)\right)$ the values on its map points which are again obtained via interpolation within adjacent planes and $s^{\pm}:=s_{i,j}(\pm\Delta\varphi)$ the associated lengths along field line. The penalization value is then prescribed as:

[TABLE]

The first terms set the actual boundary condition and use the values obtained from the field line map towards the interior domain. The second terms ensure for toroidal limiter geometries a continuous transition between both limiter sides.

A special treatment is needed for penalisation of the potentials. After having evolved the density, the electron temperature and the vorticity from time step $t$ to $t+1$ (see section 3.3) the electrostatic potential $\phi^{t+1}$ is computed at time step $t+1$ according to:

[TABLE]

yielding $\phi=\Lambda T_{e}$ in the penalisation region, where $\chi\approx 1$. The insulating sheath boundary condition implies that the generalised electromagnetic potential has to be penalised to $\psi_{\parallel}=0$, which is realised by adding the corresponding penalisation term to Ohm’s law (4). After having evolved the generalised electromagnetic potential in time to $\psi_{\parallel}^{t+1}$ the electromagnetic potential $A_{\parallel}^{t+1}$ is obtained according to:

[TABLE]

from which the parallel current is computed according to $j_{\parallel}^{t+1}=-\nabla_{\perp}^{2}A_{\parallel}^{t+1}$. Equations (15) and (16) are Helmholtz equations for $\phi^{t+1}$ and $A_{\parallel}^{t+1}$ that are solved in GRILLIX via a multigrid solver (see section 3.4). Finally, we note that we also tried slightly different methods for penalisation of the potentials [23], and the results seemed not to depend strongly on the details of the technique employed.

3.3 Time stepping

The equations are advanced in time with the 3rd order Karniadakis scheme [39]. Only the penalization term that is directly proportional to the quantity itself is treated fully implicit in time, i.e. the equations written in the form of eq. (13) are discretised in time according to:

[TABLE]

with $a_{0}=18,\,a_{1}=-9,\,a_{2}=2$ and $b_{0}=3,\,b_{1}=-3,\,b_{2}=1$. The solution for $f^{t+1}$ is trivial as the implicit penalization term on the left hand side is diagonal.

3.4 Elliptic solver

In order to compute the electrostatic potential $\phi$ from eq. (15) and the parallel electromagnetic potential $A_{\parallel}$ from eqs. (16) two Helmholtz type equations have to be solved in each time step within each poloidal plane:

[TABLE]

with given right hand side $b$, and coefficients $c_{1},\,c_{2}$ and $c_{3}$. In the global model the coefficients have in general a spatio-temporal dependency, where direct solvers become very inefficient as a costly matrix LU-decomposition would have to be performed in each time step. An efficient solution technique for eq. (17) is provided by geometric multigrid methods [40], which is implemented in GRILLIX based on a damped Jacobi smoother with trivial restriction and bilinear prolongation.

3.5 Verification

One of the main new features in GRILLIX are electromagnetic and electron inertial effects in Ohm’s law, which gives rise to shear Alfvén dynamics. The core model for the shear Alfvén wave is obtained by linearizing equations (1-6) in the isothermal limit $(T_{e}=1)$, neglect curvature, parallel ion velocity $(u_{\parallel}=0)$ and parallel resistivity $\eta_{\parallel 0}=0$. In this limit a wave equation is obtained [41] for the ‘non-adiabaticity‘ $(\tilde{n}-\tilde{\phi})$, where the tilde denotes a fluctuating quantity.

[TABLE]

with v${}_{SAW}^{2}=(1+k_{\perp}^{2})/(\beta_{0}+\mu k_{\perp}^{2})$ the phase velocity of the shear Alfvén wave, where $k_{\perp}$ is the perpendicular mode number. In the limit $k_{\perp}\ll 1$ the wave propagates at the Alfvén speed $v_{A}=\beta_{0}^{-1/2}$ and in the limit $k_{\perp}\gg 1$ at the electron thermal speed $v_{Te}=\mu^{-1/2}$. In order to verify the implementation of the electromagnetic and electron inertia effects we perform simulations with GRILLIX in a 3D periodic slab $(\mathcal{C}=0)$ without the parallel momentum equation (3) and electron thermal equation (5), but set $T_{e}=1$ and $u_{\parallel}=0$. Otherwise we run the global version of the code but initialize the density with constant background plus a fluctuation of small amplitude $\frac{\tilde{n}}{n_{bck}}=0.1$ being a mode structure in the perpendicular plane and a Gaussian along the parallel direction. The phase velocity of the divergent wave along the magnetic field is measured and compared to the analytic prediction for $v_{SAW}$. The result in fig. 2 shows an excellent agreement between GRILLIX simulations and the analytic prediction.

The presence of shear Alfvén dynamics is also numerically beneficial, as it limits the parallel electron motion resulting in a Courant-Friedrichs-Lewy number (CFL) based on the Alfén speed. In the electrostatic case $(\beta_{0}=0,\,\mu=0)$ parallel electron motion would only be hindered by parallel resistivity. A rough guess on the time step limitation in the electrostatic case can be obtained by considering the linearised vorticity equation (2) and Ohm’s law (4) neglecting electron and density variations, i.e. $n=1,\,T_{e}=1$. For a single perpendicular mode $\nabla_{\perp}^{2}\rightarrow-k_{\perp}^{2}$ a parallel diffusion equation for the electrostatic potential is finally obtained:

[TABLE]

Treating this problem explicit in time would result in a time step limitation of $\Delta t\lesssim\Delta s^{2}k_{\perp}^{2}\eta_{\parallel}$, where $\Delta s$ is the parallel grid distance. As this time step limitation is severe for large modes at low resistivity, it would make an implicit treatment practically inevitable [9, 42].

A verification of the full system via the Method of Manufactured Solution (MMS) [24] was performed in circular geometry with purely closed flux surfaces. The same setup and procedure as described in [11] for the previous simplified set of equations is applied here to the global model. The analytic MMS functions are prescribed for each dynamical field as a product of radial ($k_{\rho}$), poloidal ($k_{\theta}$ with phase shift $\delta\theta$), toroidal ($k_{\varphi}$ with phase shift $\delta\varphi$) and temporal ($\omega$) modes (see table 1). The analytic MMS functions are quite general for GRILLIX, as its numerical approach is independent of flux surfaces. The numerical error of the MMS analysis for all dynamical fields in dependence of resolution is shown in fig. 3 and follows a second order convergence, which is a good indication for correct implementation of the equations in GRILLIX.

Finally, we want to note that also a validation in slab geometry based on experiments in the Large Plasma Device (LAPD) was performed for which results can be found in [18, 19].

4 Simulation results

4.1 Setup

The parameters for our simulations are motivated from experiments with deuterium plasma in the COMPASS tokamak [20], where we normalise density and temperature to upstream separatrix values:

[TABLE]

with $a_{min}$ the minor radius. This corresponds to the following dimensionless input parameters for GRILLIX:

[TABLE]

The strongest time step limitation stems from the parallel electron heat conduction, and by linearizing this term we may anticipate a rough scaling for the time step of $\Delta t\lesssim\frac{n}{\chi_{\parallel 0}T_{e}^{5/2}}\Delta\varphi^{2}$ with $\Delta\varphi$ the toroidal grid distance between adjacent planes. In order to perform simulations at a larger time step we employ in our simulations presented here mostly a reduced heat conductivity of $\chi_{\parallel 0}=20$, but still retain the parametric dependency on temperature ($\propto T_{e}^{5/2}$). Correspondingly, we also use a reduced value for the effective sheath transmission factor of $\gamma_{e}=0.15$. In section 4.5 we investigate the effect of this by comparing to a run with more realistic heat conductivity. Development towards relaxation of the time step limitation by an implicit treatment of parallel electron heat conduction is targeted for future work.

We performed simulations in circular geometry with toroidal limiter and diverted geometry at otherwise comparable parameters. The background magnetic field for the diverted geometry is given in terms of an analytic flux function $\psi(R,Z)$ from the class of solutions described in [43] with parameters chosen as described in [44]. We define as normalised radial coordinate $\rho=\sqrt{\frac{\psi-\psi_{0}}{\psi_{X}-\psi_{0}}}$, where $\psi_{0},\,\psi_{x}$ is the poloidal magnetic flux at magnetic axis respectively at separatrix. For circular geometry we define equivalently the normalised radial coordinate as $\rho:=\frac{\sqrt{(R-R_{0})^{2}+Z^{2}}}{a_{min}}$ and the magnetic field is given in terms of a prescribed safety factor profile $q(\rho)$. The setup for the simulations in terms of radial view is shown in fig. 4 and the characteristic functions $\chi$ prescribing the location of the sheath via penalisation is illustrated for both geometries in fig. 5.

The simulations are driven via a particle source of the form

[TABLE]

and an energy source $S_{T}$ of the same form. The sources are located near the inner boundary at $\rho_{src}$ and drive the zonal averaged density $\left\langle n\right\rangle_{\rho}$ and temperature $\left\langle T_{e}\right\rangle_{\rho}$ towards prescribed values $n_{target}$ and $T_{e,target}$ within a narrow region $(w_{src})$. The parameters $c_{n},\,c_{t}$ control the rate of the sources. Outside the source region the profiles relax freely and develop self-consistently. We preferred this form for the sources to a purely flux driven source as it allows an effective control over keeping simulations within a desired parameter regime. Eventually, it is comparable to other approaches, who source their simulations via penalising the profiles near the core [23] or use a feedback control loop [8].

The main simulations analysed in section 4.2 to 4.4 were run with $32$ poloidal planes, perpendicular resolution of $h=1\,[\rho_{s0}]$ corresponding to $0.066\text{ cm}$ and a time step of $\Delta t=5\cdot 10^{-5}\,[R_{0}/c_{s0}]$. The total number of grid pints were $\approx 5.0\cdot 10^{6}$ grid points for the circular case and $\approx 14\cdot 10^{6}$ for the diverted case. A sixth order hyperviscosity $(\nu_{\perp,f}\nabla_{\perp}^{6})$ is applied in the perpendicular direction and regular diffusion $(\mu_{\parallel,f}\nabla\cdot\left(\mathbf{b}\nabla_{\parallel}\right))$ in the parallel direction, where the coefficients were chosen as $\nu_{\perp,f}=10$ and $\mu_{\parallel,f}=0.025$ cutting off turbulent spectra by smoothing structures on the grid scale. The independence of the results from these numerical parameters was checked at the circular case (see section 4.5). The radial extent of the simulations in circular geometry is $\sim 1.5\text{ cm (edge)}\,+\,3.5\text{ cm (SOL)}$ and in diverted geometry $\sim 1.5\text{ cm (edge)}\,+\,2.3\text{ cm (SOL)}$ at outboard mid plane respectively $\sim 3.2\text{ cm (edge)}\,+\,5.5\text{ cm (SOL)}$ at inboard mid plane.

The simulations were initialized with uniform background in density and temperature $(0.2)$ plus small random noise $(0.01)$. Particles and energy are injected via the sources, and the simulations enter saturated phase, which is independent of the initial state (for the circular case around $t\approx 30$) from where data is collected for performing statistical analysis. The overall simulation time for the circular case was up to $t=77$, corresponding to $\approx 1\text{ms}$. The simulations were carried out on the Marconi-A2 (KNL) partition on 16 nodes (2 MPI processes times 34 cores per node). Within 24 hours GRILLIX ran a normalised time interval of $\approx 10\,[R_{0}/c_{s0}]$ for the circular case and $\approx 4\,[R_{0}/c_{s0}]$ for the diverted case.

4.2 Circular geometry

Snapshots of density, electron temperature, electrostatic potential and parallel velocity for the circular geometry with toroidal limiter are shown in fig. 6. There is a clear difference between the closed field line region and the SOL, which is dominated by the Bohm boundary condition for the parallel velocity. Blob-like structures in the density can be observed around the last closed flux surface.

Time traces of pressure during the saturated state at low (LFS) and high (HFS) field side mid plane are shown in fig. 7. Especially the LFS-signal is highly intermittent with fluctuations of up to $200\%$, implying that a global model, which does not rely on a splitting into fluctuations and background, is indeed important. The fluctuation level on the HFS is slightly lower which results from the ballooning character of the turbulence where curvature acts stabilizing at HFS and destabilizing at LFS. It has been found that turbulence in the SOL of limited plasmas is driven by resistive ballooning modes [45, 46] with linear growth rates somewhat smaller than the interchange growth rate $\gamma_{I}=\sqrt{2R_{0}/L_{p}}$, where $L_{p}$ is the background pressure gradient length. Via the autocorrelation we may gain some insight into the characteristic time scales of the turbulent dynamics, which is computed discretely and in normalised form as:

[TABLE]

with $f(t_{n})$ the signal at discrete time point $t_{n}$. The autocorrelation function for the pressure at LFS in the region of strongest pressure gradient is shown in fig. 8, and a correlation time of $\tau_{c}\approx 0.125$ is obtained. The characteristic turbulent time scales are therefore slightly larger than the interchange time scale $t_{I}=\gamma_{I}^{-1}\approx 0.08$, where the pressure gradient length in the edge has been estimated from the self-consistently obtained profiles (see fig. 9) as $L_{p}/R_{0}\approx 0.012$. The turbulence is therefore compatible to be driven by resistive ballooning modes.

Furthermore, we analyse the simulation by computing profiles $\left\langle f\right\rangle$, fluctuation levels $\delta f=\left\langle f^{2}-\left\langle f\right\rangle^{2}\right\rangle^{1/2}$ and skewnesses $\gamma f=\left\langle f^{3}-\left\langle f\right\rangle^{3}\right\rangle/\delta f^{3}$, where angular brackets denote an average over toroidal direction and time within the saturated phase. Our statistical analysis was robust against averaging over different time windows. In fig. 9 profiles, fluctuation levels and skewnesses taken at outboard mid plane position are shown for density, electron temperature and pressure. There is a kink in the density and pressure at the last closed flux surface and a little distance outside also for the temperature. The electrostatic potential follows in the SOL $\phi\approx\Lambda T_{e}$ and deviates from this in the closed field line region where the potential at the inner limiting flux surface is prescribed as $\left.\phi\right|_{\rho_{min}}=\left\langle\Lambda T_{e}\right\rangle_{\rho=1}$. The absolute fluctuation level for the density in the closed field line region is below $0.13$ and still around $0.05$ in the near SOL. This translates into relative fluctuation levels ($\delta f/\left\langle f\right\rangle$) of around $30-40\%$, which substantiates again the importance of a global model. The skewnesses of the density and pressure approach zero in the region of its strongest gradients, which is an indication for a Gaussian probability density function and uncorrelated turbulence, which is driven in this region. Towards the core the skewnesses become negative which is an indication for holes and it becomes positive in the SOL indicating the presence of blobs.

4.3 Comparison with diverted geometry

Snapshots for the diverted case are shown in fig. 10, where again a clear distinction between closed field line region, SOL and private flux region in the dynamics is obvious. From the snapshots there seems to be a qualitative difference with respect to circular geometry: The turbulence in the saturated state is generally more quiescent especially in the SOL and we do not identify blobs at outboard mid plane as clearly as in circular geometry.

Firstly, we consider again time traces of pressure at different poloidal positions on the separatrix in fig. 11. In comparison to circular geometry (see fig. 7) the dynamics is generally more quiescent with smaller fluctuations, poloidal asymmetries are much stronger pronounced and the course into saturation is more complex and takes much longer. Until $t\approx 30$ the turbulence exhibits a strong ballooning character with violent fluctuations at the LFS, whereas the HFS is rather quiescent. After $t\approx 30$ the fluctuation levels approach each other whereas the signal at the X-point remains always rather quiescent.

Profiles, fluctuation levels and skewnesses for the diverted simulation are shown in fig, 12. Compared to the circular case (see fig. 9), the profiles are steeper and the SOL-width is smaller such that the pressure in the SOL is lower. The fluctuation levels, especially in the SOL, are also reduced with a relative level of around $15-25\%$. The skewnesses are qualitatively similar indicating that blobs in the SOL are still present, although it might not directly be visible from the snapshots in fig. 10.

To illustrate the course of the diverted simulation into its saturated state we show a series of density snapshots in the outboard mid plane region in fig. 13 at different times during the simulation. At an early time ($t=18$) -but still well after the initial onset of turbulence- there are strong fluctuations and blobs can clearly be observed. The fluctuations reduce gradually and slowly until in the saturated state the dynamics becomes relatively quiescent. This transition and the qualitative difference to the circular geometry correlates with the rise of a poloidal shear flow illustrated in fig. 14. When the flow profile of the diverted geometry at $t=18$ resembles the saturated flow profile in circular geometry the turbulence looks similar (Compare fig. 13 left with fig. 6 upper left inset). In the saturated state, a strong poloidal flow is built up exhibiting a slightly larger and qualitatively different shearing rate that swamps away turbulent fluctuations. However, we have to stress here again that that the radial electric field and therefore the flow profile in this case is probably not very realistic as the ion temperature dynamics, which has a significant impact, is not yet included in the model and the core boundary condition for the potential is also chosen somewhat ad-hoc. Besides revealing the need for an accurate and self-consistent description for the radial electric field, the study shows that geometry has a significant qualitative effect on turbulence.

Another qualitative difference between circular and diverted geometry are poloidal asymmetries. Whereas the fluctuation level in circular geometry varies within a flux surface only by a factor of two it varies in diverted geometry by an order of magnitude (see fig. 15). Approaching the X-point region from LFS there is a strong drop in the fluctuation level. The poloidal asymmetry can also be seen by comparing again the time traces in figs. 7 and 11. Especially during earlier phases of the diverted simulation the fluctuations seem to be concentrated at the LFS, whereas the HFS and especially the X-point region are rather quiescent.

To give a reason for the stronger poloidal asymmetries in diverted geometry we investigate the effect of magnetic geometry on turbulent fluctuations and consider the local magnetic shear:

[TABLE]

where $\tan\theta:=Z/(R-R_{0})$ is the geometric poloidal angle. A plot of the local magnetic shear for flux surfaces just inside the separatrix is shown in fig. 16. Whereas the local shear is obviously constant in circular geometry, it follows a complicated course in diverted geometry, i.e. it is very low in the outboard mid plane region, increases towards the top and high field side region and approaches a singularity at the X-point. Turbulent structures, which are driven in the outboard mid plane region due to unfavourable curvature, become distorted in the perpendicular direction due to local magnetic shear (see also [47]). Being strongly distorted, i.e. especially in the vicinity of the X-point, the fluctuations are damped due to perpendicular dissipation. Therefore the X-point ultimately acts as kind of barrier for fluctuations [21] (See also resistive X-point mode by Myra et al. [48]). This explains the drop of the fluctuation level near the X-point towards HFS in our simulation (see fig. 15b). There is also a drop in fluctuation level in the top region towards HFS, where local magnetic shear is also relatively large.

In contrast to flux-aligned approaches GRILLIX does not suffer from coordinate singularity and a loss of resolution near the X-point due to flux expansion, but the FCI approach allows an accurate treatment of the dynamics around the X-point. A detailed view of the pressure at the X-point is shown in fig. 17. Fluctuations approaching the X-point fan out radially becoming ever narrower in the poloidal direction, which illustrates the mechanism described in the previous paragraph.

4.4 Comparison to Boussinesq system

Although the Boussinesq approximation is factually not justified in the edge and SOL, it has often been employed in various codes and in various forms for numerical reasons (e.g. [11, 7, 8, 9]). For that matter the effect of the Boussinesq approximation has been a long standing discussion, and since recently there have been several developments on its abolishment [6, 10]. At the same time the impact of the Boussinesq approximation has mainly been studied at isolated phenomena, i.e. blobs [6, 16, 17], whereas we investigate here the effect of the Boussinesq approximation in a fully turbulent environment. For this purpose we modify the vorticity equation (2) to:

[TABLE]

where we define the Boussinesq vorticity as $\Omega_{bsq}=\nabla\cdot\left(\frac{1}{B^{2}}\nabla_{\perp}\phi\right)$, i.e. we drop the spatio-temporal dependence of the density in the polarisation term completely. We note that we also tried another form for the Boussinesq approximation, where the density is taken out of the divergence with retention of its spatio-temporal dependence, i.e. $\Omega_{bsq(alt)}=n\nabla\cdot\left(\frac{1}{B^{2}}\nabla_{\perp}\phi\right)$. This alternative form breaks the conservation property of the quasi-neutrality equation, for which reason we observed a strong spurious poloidal flow rising that eventually caused our simulations to crash [18, 19].

Based on the circular limited case, we directly compare our simulation results from section 4.2, obtained with the global model, to the Boussinesq reduced model. We restarted the circular case from scratch, and except from the Boussinesq modification there are no other modifications neither in physical nor numerical parameters. The result of this comparison is shown in fig. 18, i.e. profiles, fluctuation levels and skewnesses at outboard mid plane. Our results were again checked to be robust against statistical noise via averaging over distinct time windows within the saturated phase. There are only small quantitative differences between the global and the Boussinesq model. There is basically no difference in the temperature profile and only a small difference of around $10\%$ in the density profile around the separatrix. The fluctuation levels differ at most around $10\%$ and the skewnesses are very similar in both cases. In conclusion, the Boussinesq approximation has only minor quantitative effects for the case considered here.

4.5 Convergence analysis and impact of heat conductivity

A posteriori, we subject our results to a convergence check and study the impact of heat conductivity, as we decreased it artificially. We perform this study only for circular geometry as it is computationally generally less expensive than the diverted case and its saturation time is shorter.

For a convergence check we ran the circular case at nominal resolution, at half resolution and at a resolution that was increased by a factor of one third. We correspondingly adapted also other numerical parameters, i.e. we decreased numerical dissipation coefficients with increasing resolution. Moreover, we also varied the poloidal and radial decay lengths of the penalisation functions $\chi$ in order to exclude spurious effects from our penalisation approach. The numerical parameters employed are listed explicitly in table 2. The obtained profiles and pressure fluctuation level at LFS are shown in fig. 19. Whereas there is a deviation in the temperature profile to the coarse resolution case, the fine and nominal resolution match here very well. The density profiles between nominal and fine resolution match overall well apart from the density shoulder at the last closed flux surface, which is slightly more pronounced at finer resolution. Also the pressure fluctuation converges and deviation between nominal and fine resolution could also be owed to the fact that for the fine case not as much statistics as for the other cases was available due to computational constraints.

Compared to realistic COMPASS parameters we employed a significantly reduced artificial heat conductivity of $\chi_{\parallel 0}=20$ instead of $\chi_{\parallel 0}=340$ and correspondingly also decreased the effective sheath transmission factor to $\gamma_{e}=0.15$ instead of $\gamma_{e}=2.5$. In order to study the effect of this we restarted with reduced time step the circular limited simulation from a saturated state with a more realistic heat conductivity of $\chi_{\parallel 0}=140$ and $\gamma_{e}=1.0$. The obtained profiles and pressure fluctuation level are shown in fig. 20. While the increased heat conductivity does not alter the density profile significantly it steepens the temperature profile in the SOL. The fluctuation level reduces slightly which is consistent with [46], where a similar study was carried out.

5 Conclusions and Outlook

By incorporating full parametric dependencies and relaxing the Boussinesq approximation, GRILLIX was extended to a global fluid turbulence code for the tokamak edge and SOL, i.e. no assumption about fluctuation amplitudes of density or temperature is made. Further new features are electromagnetic and electron-thermal dynamics, and the implementation of the extended model was verified by analytical means and the Method of Manufactured Solutions (MMS). The flux-coordinate independent approach (FCI) is employed in GRILLIX in order to deal with realistic geometries avoiding coordinate singularities at the X-point or separatrix. As boundary contours do in general not conform with the computational grid nor the magnetic field as exceptional direction, a penalisation technique is used in order to treat boundary conditions at the target plates.

Turbulence simulations in circular geometry with toroidal limiter and in realistic diverted geometry at otherwise comparable parameters were presented, where parameters were chosen being characteristic for the COMPASS tokamak. A relative fluctuation level of around $30-40\%$ in the SOL, with isolated highly intermittent phenomena of up to $200\%$ fluctuation level were observed in circular geometry, which implies that a global description for SOL turbulence is indeed important. At the same time a direct comparison to a Boussinesq reduced model revealed that the Boussinesq approximation had only minor quantitative effects. However, this conclusion might possibly not hold true for other parameter regimes [19] and therefore it is safest to abolish it consequently. Besides this, other global features, e.g. parametric dependencies of the heat conductivities and resistivity, might play a more important role. In diverted geometry the turbulence was found to be generally more quiescent in the saturated phase, which is related to a different zonal flow structure. Moreover, the diverted geometry exhibits stronger poloidal asymmetries that can be attributed to local magnetic shear, which is fundamentally different in the edge region between circular and diverted geometry. Approaching regions of strong local magnetic shear turbulent structures become distorted towards the high field side and subsequently subject to enhanced perpendicular dissipation. As this mechanism is very strong near the X-point, it tends to disconnect low field side and high field side, where curvature is favourable [21]. In conclusion the presented studies point out the importance of global effects due to the presence of large intermittent fluctuations and qualitative differences between circular and diverted geometry.

However, the results have to be considered qualitative as there are yet important effects in the physical modelling missing, i.e. ion thermal dynamics, interaction with neutrals and more realistic sheath boundary conditions that could also take into account glancing angles of incidence. Our study revealed also the importance of the radial electric field, and its self-consistent description requires either an adjustment of the core boundary conditions or simply a complete abolishment of the core boundary by performing computationally intense full tokamak simulations. Besides this, there are also numerical constraints in order to achieve realistic parameter regimes concerning computational efficiency. The strongest time step limitation arises from the parallel non-linear heat conduction scaling strongly with $\chi_{\parallel}\propto T_{e}^{5/2}$, for which an implicit treatment will be necessary. Efforts to extend GRILLIX in these directions are currently ongoing.

In the near future our studies would concentrate on investigation of advanced divertor concepts such as double-null, snowflake or super-X configurations, whose treatment is straight forward with GRILLIX as its numerics is independent of flux surfaces.

6 Acknowledgements

The authors would like to thank A. Bottino for computational support. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.