Landau-fluid simulations of edge-SOL turbulence with GRILLIX Open Access

Fluid models are still the workhorse for plasma edge turbulence simulations of magnetic confinement fusion devices. Codes, such as BOUT++,¹ GDB,² GBS,³ SOLEDGE3X,⁴ or GRILLIX,⁵ are designed to solve fluid equations in the edge and scrape-off layer (SOL) of magnetic fusion experiments, with their challenging geometry. Any fluid model requires a closure to become applicable. The most common choice is the Braginskii closure⁶ that is rigorously valid in the short-mean-free-path regime, i.e., in a highly collisional plasma. Conditions violating this assumption can already be present in the near SOL of present-day experiments, where the Spitzer–Härm formula,⁷ used in the Braginskii closure, can result in vast overestimation for the parallel heat conductivity.

To study reactor relevant conditions, it is, therefore, necessary to extend the model into the direction of low collisionality. A common way to deal with the overestimation of the parallel heat conductivity by the Spitzer–Härm formula is the usage of heat-flux limiters.^8–10 Those limiters have the disadvantage of introducing a free parameter into the model that has to be tuned and cannot be estimated easily in advance. Furthermore, the radial turbulent transport depends strongly on the parallel heat conductivity and, therefore, on this free parameter.^8,11

A more elaborate approach to predict parallel conductive heat fluxes for low collisional plasmas is the Landau-fluid closure.¹² Hammett et al.¹³ proposed this closure for collisionless plasmas by introducing linear Landau-damping¹⁴ into fluid models. It was extended for arbitrary collisionality by Umansky et al.¹² as suggested by Snyder et al.¹⁵ and Beer et al.¹⁶ The Landau-fluid closure was originally formulated in wave number space, but most edge-SOL turbulence codes act in configuration space due to geometrical constraints. As a solution, an efficient fast non-Fourier method was presented in Dimits et al.¹⁷ 

We illustrate the Landau-fluid closure within a 1D setup¹⁸ and compare it with the Braginskii closure⁶ and the collisionless Hammett–Perkins closure.¹³ All models are tested in various regimes of collisionality and significant differences in amplitude and in the form of the heat fluxes are found. The overestimated heat conductivity predicted by the Braginskii closure is reproduced and for the other two closures non-local heat fluxes are observed.

The Landau-fluid closure is implemented into the edge plasma turbulence code GRILLIX,⁵ which employs the flux-coordinate-independent (FCI) approach.^19,20 Calculating the heat flux predicted by the Landau-fluid closure requires solving a set of decoupled elliptic equations along magnetic field lines. In globally field-aligned approaches, this becomes a 1D problem that can be solved efficiently, as, e.g., shown with BOUT++.²¹ However, within the FCI approach, this is a fully 3D problem, which is computationally more demanding. The implementation in GRILLIX is verified via the method of manufactured solutions.²² 

For the first time, a non-linear turbulence simulation in realistic divertor geometry with the Landau-fluid closure is presented and compared with simulations using the Braginskii closure. Three simulations are set up, the first with Braginskii closure and a strong heat-flux limiter, the second one also with Braginskii closure but a weak heat-flux limiter, and a third one employing the Landau-fluid closure. All three simulations also use recently implemented magnetic flutter effects and neoclassical ion viscosity. Clear differences between the simulations are found, and most pronounced is an ion temperature asymmetry for the strongly limited Braginskii case due to heat-flux limiters. A detailed investigation of the Landau-fluid closure during a turbulence simulation is performed and the computational performance of the model is assessed. We show that non-local heat fluxes are present and relevant for turbulence simulations. Therefore, the Landau-fluid closure cannot be mimicked with a flux-limited closure, regardless of the choice of the free flux-limiting parameter.

This paper is organized as follows. In Sec. II, the necessary theory is explained, Sec. III deals with the implementation of the Landau-fluid closure in a 1D setup and in the FCI framework of GRILLIX. In Sec. IV, three turbulence simulations are presented and analyzed. Finally, in Sec. V, a summary and an outlook are provided.

Fluid moment equations are generally obtained by taking velocity moments of the Boltzmann or Vlasov equation. In this infinite set of fluid equations, the full information of the initial kinetic system is still preserved. However, to make a fluid model applicable, a closure approximation is needed to truncate the fluid hierarchy. Since the truncation of the equation system is already an assumption, there is no single correct solution to close the system. Therefore, different choices are available. The set of fluid equations used in GRILLIX can be found in  Appendix A. For this set of equations, we need a closure for the parallel conductive heat flux $q ∥$ and for the stress tensor $Π$⁠. The closure for $Π$ is briefly discussed in  Appendix A. The focus of this work is set on the closure for $q ∥$ and from now on the term fluid closure always refers to a closure for the parallel conductive heat flux.

For decreasing collisionality, however, the closure loses its validity. In regimes of high collisionality, collisions can drive the plasma distribution function back to a Maxwellian very efficiently; but, this mechanism weakens as the plasma becomes less collisional. Specifically, the issue that arises in the equation for the closure is the dependency $κ ∥ j BR ∝ T j 5 / 2$⁠. Only a very collisional plasma obeys this temperature dependence; however, for low collisional plasmas, the heat conductivity is dramatically overestimated with this closure.

We calculate the magnitude of the peak heat fluxes for T₀ = 3, 30, and 300 eV and plot them in Fig. 3. The three heat flux models align with our expectations by comparing Fig. 3 with Fig. 1. For this setup, T₀ = 3 eV is a highly collisional case, T₀ = 30 eV an intermediate case, and T₀ = 300 eV is nearly collisionless. To investigate the Landau-fluid closure in more detail, we take a closer look at the case with T₀ = 30 eV. In this regime of intermediate collisionality, we examine the right half of our domain, i.e., $x ∈ [ 0 , 0.5 ]$⁠. In Fig. 4, we show the initial temperature profile, the corresponding Braginskii, Hammett–Perkins, and Landau-fluid heat flux and the contribution of every single Lorentzian to the Landau-fluid heat flux. Here, we confirm the observation of non-local behavior of the latter two heat fluxes as already reported in Chen et al.¹⁸ By non-local we mean heat fluxes in regions where no local temperature gradient is present and a local heat flux model as the Braginskii closure would predict no heat flux, in Fig. 4, this corresponds to x > 0.3.

For the single Lorentzians, we observe that for this case mainly numbers 2, 3, and 4 are important. Numerous tests with different temperatures and initial temperature distributions revealed two general trends. First, as collisionality decreases, higher numbered Lorentzians become more important until the highest number is reached and, therefore, the limit of the model. One can think of the bright spot moving upwards in the plot and eventually out of the limits. When this limit is reached, the model starts to underestimate the heat flux strongly as there are no Lorentzians left to represent it. For this simple test case, the model with seven Lorentzians performs well up to $T 0 = 1 keV$⁠. This value is above typical temperature values in the SOL, but 1 keV can be exceeded in the outer closed-field-line region of present-day devices. Second, for this simple 1D setup always just a few neighboring Lorentz functions are important, as can be seen in Fig. 4 as well, and out of them the lower numbered Lorentzians are responsible for the non-local behavior.

The equations that need to be adapted in GRILLIX are the two temperature equations [Eqs. (A5) and (A6)], in which the Braginskii expressions [Eq. (1)] for $q ∥ e$ and $q ∥ i$ are simply replaced by the Landau-fluid expressions [Eq. (7)]. In a globally field-aligned approach, Eq. (7) would simply reduce to a set of independent 1D problems, if only the unperturbed magnetic field $B 0$ was taken into account for parallel derivatives.²¹ However, in the FCI^20,28 approach employed by GRILLIX, the set of field-aligned elliptic equations (7) are 3D problems since there is no global field-aligned coordinate in FCI. These 3D problems are solved using the PIM (parallel iterative methods) library,²⁹ which provides an efficient GMRES (Generalized Minimal RESidual) algorithm to solve the global problem that is distributed over different MPI processes with matrix-free methods. To speed up convergence, a Jacobi preconditioner is used.

Boundary conditions in GRILLIX are divided up into perpendicular and parallel. The perpendicular ones are homogeneous Neumann boundary conditions for parallel heat flux. No major influence of this boundary condition is expected and found since buffer zones, particle, and heat sources are active in those regions as well. The parallel boundaries are enforced as source terms via an immersed boundary method.^5,30 A penalization function χ is introduced that has a value of χ = 0 inside the domain, χ = 1, where the boundary conditions are applied with a penalization strength $1 / ε$ and a smooth transition between the two domains. With this penalization method, we have to modify Eq. (7) and write down the boundary conditions directly into the corresponding equations.

The parallel boundary conditions for heat flux at the divertor targets are sheath boundary conditions, which are stated in Eq. (13). For the Braginskii closures, this boundary condition has to be applied as a Neumann boundary condition on the temperature. On the other hand, with the Landau-fluid closure, we can apply the sheath boundary condition directly to the heat flux as a Dirichlet boundary condition.

In GRILLIX, the magnetic field is normalized by B₀, the density by n₀, the temperatures via $T e 0$ and $T i 0$⁠, the electrostatic potential by $T e 0 / e$⁠, the parallel velocities by sound speed $c s 0 = T e 0 / M i$⁠, the current by $c s 0 e n 0$⁠, the perpendicular length scales by $ρ s 0 = c T e 0 M i / ( e B 0 )$⁠, the parallel length scales by R₀, the time by $R 0 / c s 0$⁠, and the parallel component of the magnetic vector potential by $β 0 B 0 ρ s 0$⁠, with $β 0 = 4 π n 0 T e 0 / B 0 2$⁠.

We investigate and compare the Landau-fluid closure with the previously used Braginskii closure in a fully non-linear turbulence simulation resembling the ASDEX Upgrade experiment. Three simulations are performed, the first one using the Braginskii closure with a strong heat-flux limiter (⁠ $α = 0.1$⁠), the second one also with the Braginskii closure but a weaker heat flux limiter (⁠ $α = 1.0$⁠), and the third one employs the Landau-fluid closure. The complete physical model is given in  Appendix A, for a comprehensive discussion of the model, we refer to Zholobenko et al.³¹ In addition to the Landau-fluid closure, two further model improvements have been applied since then, namely, magnetic flutter and neoclassical ion viscosity. Magnetic flutter, as introduced in Eq. (9), is used for all parallel operators in our system and for the Landau-fluid equations (10) and (11) just in the temperature gradient on the right-hand side, after preliminary tests showed that including it in the left-hand side does not change the physical results, but the computational cost of the 3D-solver is increased. The neoclassical ion viscosity extension is briefly discussed in  Appendix A. Both extensions will be discussed in detail in separate publications. The geometry and the equilibrium of these simulations are based on ASDEX Upgrade discharge #36190, an attached L-mode with a plasma current of 800 kA, q₉₅ = 4.4, and an average triangularity of $δ = 0.21$⁠. The toroidal magnetic field strength is $B tor = 2.5 T$ on axis in favorable configuration, i.e., the ion- $∇ B$ drift points toward the active X-point. The plasma was heated in the experiment with Ohmic heating and neutral beam injection. After subtracting the radiation losses, which are not modeled by GRILLIX, the total input power was 475 kW in the experiment. Two species are considered in the simulation, electrons, and deuterium ions. In all simulations, the density, electron, and ion temperature are kept constant at the inner core boundary of our simulation domain by an adaptive source at $ρ pol = 0.91$⁠, which maintains the values $n ( ρ pol = 0.91 ) = 2.0 × 10 19 m − 3$ and $T e ( ρ pol = 0.91 ) = T i ( ρ pol = 0.91 ) = 350 eV$⁠. Additionally, a diffusive neutrals model is used for all the simulations,¹¹ which requires a fixed neutrals density at the divertor as a boundary condition, this value is here chosen to be $5.0 × 10 17 m − 3$⁠. All parameters remain the same for the three simulations, except for the different heat flux closures employed.

The time traces of flux-surface-averaged density, electron, and ion temperature at $ρ pol = 0.998$ for all three simulations are depicted in Fig. 5. Both the Braginskii case with weak limiter (⁠ $α = 1.0$⁠) and the Landau-fluid case seem to approach convergence. However, the Braginskii case with strong limiter (⁠ $α = 0.1$⁠) shows a second rise in the electron and ion temperature starting at t = 2.5 ms. This originates in a strong asymmetry in the ion temperature, which will be explained in more detail in Sec. IV C.

To compare the three simulations directly, we plotted the outboard-midplane (OMP) profiles of density, electron, and ion temperature on top of each other for the three simulations at t = 2.5 ms in Fig. 6, averaged toroidally and over $Δ t = 100 μ s$⁠. The density profiles exhibit distinct “shoulders” near the separatrix, especially for the two cases with Braginskii closure. These shoulders are caused by neutral gas that is ionized near the X-point. For the Landau-fluid closure, this effect is less pronounced. Since the neutrals density at the divertor is a free parameter in the model, a careful scan would be necessary to match the actual experimental values. However, this work focuses on the investigation of the different heat flux closures and their respective influence, an extensive validation against the experiment we leave for future work. The electron [Fig. 6(b)] and ion temperature profiles [Fig. 6(c)] are more directly influenced by the respective heat-flux closure. Here, we see for the temperature of both species that the Landau-fluid case falls in between the two Braginskii cases in terms of separatrix temperature. Deep in the edge region (⁠ $ρ pol < 0.97$⁠) the three cases predict similar temperature profiles, except for the ion temperature profile of the Braginskii case with $α = 0.1$⁠, where the different shape of the ion temperature profile is due to an asymmetry in ion temperature, which is examined in Sec. IV C. Looking closer especially at the temperature profiles, a decrease in the fluctuation amplitudes by going from $α = 0.1$ to $α = 1.0$ to the Landau-fluid model is visible. This indicates a change in the underlying turbulent dynamics. To investigate this a little deeper, the radial electric field is shown in Fig. 6(d). The origin and properties of the radial electric field in simulations with GRILLIX were analyzed in great detail by Zholobenko et al.³¹ for the same setup and compared with the experimentally measured values. We find that with the Landau-fluid closure the quantitative values of the radial electric field are significantly closer to the experiment.¹¹ An extensive discussion of the influence of the closure on the radial electric field will be left for the future since it is beyond the scope of this work.

Comparing the input power in Table I, we observe the Braginskii case (⁠ $α = 1.0$⁠) and Landau-fluid being similar in terms of total input power and close to the experiment. The case with Braginskii closure (⁠ $α = 0.1$⁠) deviates strongly in electron input power, regardless of the just slightly flatter electron temperature profile [Fig. 6(b)]. Nevertheless, the flux limiters work as expected since an increase in α leads to smaller temperature fluctuations, therefore, to decreased radial heat transport, so less input power is needed to keep the temperature constant at the core boundary. In Zholobenko et al.,¹¹ a very similar case was considered in the original Braginskii limit $( α → ∞ )$⁠. Such simulations are computationally challenging as the parallel heat conductivity is becoming very stiff. A further steepening of the profiles and a very low input power, much lower than typical experimental values, were obtained.

The already mentioned asymmetry in ion temperature for the Braginskii case (⁠ $α = 0.1$⁠) becomes apparent by looking at the corresponding ion temperature across the whole simulation domain in Fig. 7 in comparison with the two other cases. A strong up-down asymmetry is present in the ion temperature which is shifted in anti-clockwise direction due to the poloidal E × B rotation. Since this asymmetry is still present at the OMP, it determines the shape of the ion temperature profile [Fig. 6(c)].

The limitation of the parallel conductive ion heat flux for the Braginskii (⁠ $α = 0.1$⁠) case is indeed apparent from Fig. 9(a), where $⟨ | ⟨ q ∥ i ⟩ φ | ⟩ θ , t$ is plotted. On the other hand, the electron heat flux in Fig. 9(b) shows no big difference between the two Braginskii cases, but for the Landau-fluid case. First, we note that the limiting free-streaming heat conductivity is about a factor $M i / m e ≈ 60$ higher for electrons than ions. Second, we have to keep in mind that the heat conductivity is limited and not the heat flux directly. Therefore, we additionally plot $⟨ | ⟨ ∇ ∥ T j ⟩ φ | ⟩ θ , t$ in Figs. 9(c) and 9(d). As expected, we find for electrons in the case Braginskii (⁠ $α = 0.1$⁠) that the similar amplitude of the heat flux in comparison with the case Braginskii (⁠ $α = 1.0$⁠) can just be obtained by a significantly larger parallel temperature gradient in the quasi-stationary state. We want to highlight that the Braginskii case (⁠ $α = 1.0$⁠) is still limited in comparison to an unlimited Braginskii case (⁠ $α → ∞$⁠), which would provide significantly larger parallel heat fluxes for comparable temperature gradients. Interestingly, we note that the parallel electron heat flux for the Landau-fluid case is large despite the small parallel temperature gradient, which illustrates that no simple linear connection between $∇ ∥ T$ and $q ∥$ exists for this case. This finding is investigated in more detail in Sec. IV E.

To ensure that the asymmetry in ion temperature is caused by the curvature term and balanced by the parallel-heat-flux term, we examine those terms along the flux surface at $ρ = 0.965$⁠. The value of ρ is chosen to match the radial position of the largest temperature asymmetry. This analysis is performed for all three cases and is depicted in Fig. 8, where the black line shows the sum of all remaining terms of Eq. (A6). The ion temperature asymmetry is strongest for the case with Braginskii closure (⁠ $α = 0.1$⁠) and weakest for the case with the Landau-fluid closure. The analysis confirms that the curvature term is building up the asymmetry as it peaks for all cases near the maximum of the temperature asymmetry and is again strongest for the case with Braginskii closure (⁠ $α = 0.1$⁠). We cannot compare the three cases one to one, since every case is a full non-linear turbulence simulation and changing the parallel heat conduction changes the state of the whole system. Nevertheless, the parallel-heat-flux term in Fig. 8 for Braginskii (⁠ $α = 1.0$⁠) and Landau-fluid together with the sum of the remaining terms of the ion temperature equation (black line) counteracts the curvature term and decreases the ion temperature asymmetry. For Braginskii (⁠ $α = 0.1$⁠), we see the heat flux term being close to zero throughout the whole flux surface because of the strong heat-flux limiter. Comparing the two terms to the sum of all remaining terms of the ion temperature equation shows that the two terms are not just small perturbations, but significant contributions.

A similar asymmetry in ion temperature was observed and described already in circular geometry by Zhu et al.³³ with the GDB code. In their simulations, a weak up-down asymmetry was observed in all quantities but the most prominent one in ion temperature. By turning off the parallel heat conduction the asymmetry got more pronounced, which can be understood as an extreme case in terms of heat-flux limiters by setting α = 0.

In the theory section, we have discussed already that the Landau-fluid closure employs higher numbered Lorentzians to represent the parallel heat flux as collisionality decreases. This correlation was found with a simple 1D model. To investigate if this characteristic is also visible in a turbulence simulation, where temperature distributions along field lines are far more complex than a simple Gaussian, we plot the flux surface averaged absolute values of the individual Lorentzians $⟨ | ⟨ q ∥ j , n LF ⟩ φ | ⟩ θ , t$ over $ρ pol$ in Fig. 10. For electrons [Fig. 10(a)], it is clearly observed how the number of the relevant Lorentzian is increasing with decrease in $ρ pol$⁠. This is perfectly in line with our expectation, since with decreasing $ρ pol$ the temperature increases and therefore the collisionality reduces. For ions, the same behavior is visible although not as distinctly separated as for electrons. This method offers also a straightforward approach for testing if the number of Lorentzians used in the simulation is sufficient. If the ensemble of active Lorentzians would reach the highest number or even move out of the frame, then the number of Lorentzians is insufficient for resolving the range of collisionality. Furthermore, two additional simulations were conducted with 3 and 12 Lorentzians to ensure the consistency of this heat-flux closure. For three Lorentzians, we found that the simulation is indeed under-resolved in terms of Lorentzians. In this case, a similar ion temperature asymmetry as in the Braginskii case (⁠ $α = 0.1$⁠) was observed, since the parallel ion heat flux is artificially damped (see Fig. 11), as there are no higher numbered Lorentzians to represent it.

This mechanism becomes evident from Fig. 2, where the parallel heat flux predicted by the Landau-fluid closure approaches zero for very low collisionality. The simulation with 12 Lorentzians does not show any major differences compared with the standard case with 7 Lorentz functions.

We want to investigate whether the non-local heat fluxes predicted for Landau-fluid closure in the 1D setup (see Fig. 4) are also relevant in a 3D turbulence simulation. Furthermore, we want to investigate if the effect of the Landau-fluid closure could be approximated for this case by an appropriate local closure.

Furthermore, we exclude the area near the core boundary since the temperature and density sources are active there. Second, $α LF$ can reach high values when $κ ̃ ∥ j LF$ itself becomes very large. This is the case we are interested in since it implies a non-locality as we observe a heat flux in the region of vanishing temperature gradients.

The calculated effective $α LF$ for electrons and ions is depicted in Fig. 12. We have cut the color scale since $α LF$ is reaching excessively high values in this region. The value of the cut is arbitrary and chosen here for visualization reasons, but we want to denote that the values in the domains reach higher values (⁠ $α j LF > 1000$⁠). For both electrons and ions, we find values for $α LF$ that are far beyond the usually used limit of α = 3. This is a strong evidence that the non-local behavior of the Landau-fluid closure is not just academically relevant in simple 1D cases (Sec. III A), but also in non-linear edge turbulence simulations. The high values of $α LF$ also indicate that the parallel temperature gradients are significantly smaller with the Landau-fluid closure than with the Braginskii closures, considering that the heat fluxes are of the same order of magnitude (see Fig. 9). Furthermore, it is interesting to note that $α LF$ also reaches negative values, which implies that heat fluxes in those regions flow upward the temperature gradient. This effect is, by definition, not possible for a local closure. The conclusion of constant flux-limiters in space and time being unable to represent non-local heat-flux models was drawn in Omotani et al.³⁴ (Fig. 4) for a 1D time-dependent system already. However, we do not want to make estimations for relevant values of α to conduct simulations with a flux-limited closure, since this rather pragmatic comparison between local and non-local closures is not expected to work quantitatively. Detailed quantitative investigations of non-local heat-flux models in simpler geometries have been already conducted by Brodrick et al.³⁵ We tried to investigate the Landau-fluid closure specifically for turbulence simulations.

Another way to visualize this behavior of the Landau-fluid closure is by plotting the effective heat conductivity [Eq. (16)] for every point in the domain in Fig. 12 over the corresponding temperature in Fig. 13. We can draw the same conclusion here, since the single points spread over an area they cannot be described by an analytic regression function dependent on temperature as it is possible for the flux-limited expressions, especially for electrons.

To give an impression on the computational cost, for the Landau-fluid simulation 3.8 × 10⁶ time steps were performed for simulating well into the saturated state until 3.5 ms of physical time. This simulation was carried out on the MARCONI SKL partition from CINECA, employing eight compute nodes with 48 cores each and took approximately 40 days, i.e., 3.7 × 10⁵ CPUh in total.

As pointed out in Stegmeir et al.,³⁰ parallel dynamics play a crucial role for code performance. From the computational perspective, the Landau-fluid and Braginskii model are very similar, as the implicit treatment of the Braginskii heat-flux term also results in an elliptic problem along magnetic field lines, similar to Eqs. (10) and (11), which is solved with the same GMRES method from the PIM²⁹ library. The difference being that for the Braginskii model only one such elliptic problem has to be solved per species, whereas for the Landau-fluid model, N problems depend on the number of Lorentzians. Otherwise, we were able to use similar time step sizes for all simulations.

For a direct performance comparison, we continued all five simulations for a few thousand time steps starting from the last stored point in time on eight nodes of the RAVEN cluster, which employs Intel Xeon IceLake-SP 8360Y processors with 72 cores per node, from the Max Planck Computing and Data Facility. In Table II, the average duration for a single time step of the plasma model is given, which comprises the evaluation of explicit terms (RHS), 2D elliptic solves for the electrostatic and parallel perturbed electromagnetic potential, and 3D field-aligned solves (where the necessary MPI communication is included) for the parallel heat fluxes and parallel momentum dissipation. As expected, the Landau-fluid model comes with an expense in computational time increasing with the number of Lorentzians. The increase of 46% for 7 Lorentzians can be considered acceptable given the gained validity of the physical model and the removal of free parameters.

Finally, the bottleneck of the code shifts from the 2D solvers to the 3D solvers. While much optimization effort has already gone into the 2D elliptic solvers, for which even a GPU porting is currently ongoing, little work has been invested into the 3D solver so far, but only a naive GMRES algorithm with Jacobi preconditioning is employed. Therefore, a large span of optimization is still possible, by employing more advanced algorithms, e.g., algebraic multigrid methods offered by the HYPRE³⁶ library that comes already with GPU support.

The Landau-fluid closure for the parallel heat fluxes was implemented into the edge turbulence code GRILLIX. As illustrated in a 1D setup, this introduces non-local effects and self-consistently limits heat-flux levels to more realistic levels toward collisionless regimes, where the Spitzer–Härm formula would result in a vast overestimation. For the implementation within the FCI approach, a set of 3D elliptic problems has to be solved, whereby the incorporation of electromagnetic flutter terms is straightforward. The same solver that was already employed for the implicit treatment of the Braginskii heat conduction term could be reused. Verification was performed via the Method of Manufactured solutions.

Three turbulence simulations, two with the Braginskii closure employing different heat flux limiters and one with the Landau-fluid closure, were performed. An emerging up-down asymmetry in the ion temperature was observed, which was most prominent for the Braginskii case with strong limiters (⁠ $α = 0.1$⁠). We provide evidence that the asymmetry is caused by the strong limitation of the parallel heat conductivity that is not sufficient anymore to balance the diamagnetic heat flux driving the asymmetry. The non-locality of the Landau-fluid closure during a turbulence simulation is examined by calculating an effective flux-limiting parameter $α LF$ in the closed-flied-line region of our simulation. The excessively high values of $α LF$ indicate the presence of non-local dynamics. Moreover, the large span in the effective $α LF$ suggests that the Landau-fluid closure can hardly be mimicked with a simple flux-limiter in the Braginskii model. For the Landau-fluid simulation, we generally observe that toward lower collisionality the heat flux calculated is carried by higher-numbered Lorentzians. A convergence check in the number of Lorentz functions showed that seven Lorentzians were sufficient for the considered case. In the simulation with only three Lorentzians, the heat flux is underestimated at low collisionality, which resulted again in a strong up-down asymmetry of the ion temperature.

Concerning performance, we measure an increase of 46% in computational time of the Landau-fluid closure in comparison with the Braginskii closure for the given case. This appears acceptable considering the increased validity of the physical model and the elimination of free parameters. Moreover, there is still high potential for optimizing the 3D solver.

A comparison between the Landau-fluid closure, which includes linear Landau-damping into the fluid model, and a gyro-kinetic code as, e.g., GENE-X³⁷ would provide further insights but is left here for future work.

The simulation of this L-mode discharge shows already influences of the Landau-fluid model, even at relatively low temperatures and therefore still quite high collisionality. Still there is a good agreement in profiles between the Landau-fluid closure and the Braginskii model with appropriately chosen heat flux limiters. The importance of the new closure might increase for simulations of reactor-relevant regimes, like H-mode³⁸ or I-mode,³⁹ with significantly lower collisionality in the edge region. These operational scenarios are especially interesting to study in the edge region, as transport barriers arise there.⁴⁰ The Landau-fluid closure, which extends the model into the direction of low collisionality, will play a crucial role in these attempts and provides a step toward the goal of being able to simulate reactor-relevant scenarios. This process shall be accompanied by detailed validation of all the new features in GRILLIX against experimental data.

The authors would like to thank B. Zhu for fruitful discussions as well as K. Eder, J. Pfennig, and P. Ulbl for discussions on physical and design aspects. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The authors have no conflicts to disclose.

Christoph Pitzal: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Methodology (equal); Software (equal); Visualization (equal); Writing – original draft (lead). Andreas Stegmeir: Conceptualization (supporting); Formal analysis (supporting); Methodology (equal); Software (equal); Supervision (equal); Visualization (supporting); Writing – original draft (supporting). Wladimir Zholobenko: Conceptualization (supporting); Formal analysis (supporting); Methodology (equal); Software (supporting); Supervision (supporting); Visualization (supporting); Writing – original draft (supporting). Kaiyu Zhang: Investigation (supporting); Software (supporting); Visualization (supporting). Frank Jenko: Conceptualization (supporting); Funding acquisition (lead); Methodology (supporting); Supervision (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The neoclassical extension of the ion viscosity G essentially adapts the Braginskii coefficient $η i 0$ to $η ̃ i 0$ according to Eq. (39) in Rozhansky et al.⁴¹ This follows the work of Hirshman and Sigmar from 1981.⁴²  Fig. 1 therein shows, in particular, that this introduces an upper limit for the ion viscosity coefficient, while the Braginskii expression diverges as $T i 5 / 2$⁠, similarly to the (ion) heat conductivity. Neoclassical heat viscosity is not yet included. For the work here, we approximate the connection length as $L c ≈ q 95 R 0$⁠, and the inverse aspect ratio as $ε ≈ 0.3$⁠, these are typical values for ASDEX Upgrade. Details on the physical effect of this extension will be published separately. We have to make the remark that all simulations presented in this paper contain a small mistake in the form of the neoclassical extension. The collisionality was implemented with a temperature dependency of $ν i ∝ T i 5 / 2$ instead of the correct dependency of $T i 2$⁠. Although small quantitative changes are possible, we do not expect a qualitative change in the results, especially as the same model (including the mistake) was used for all simulations.

With these analytical solutions, source terms are calculated as described as in Salari et al.²² Those source terms are added on the right-hand side of all equations in  Appendix A and also for the two Landau-fluid equations (10) and (11). Those source terms are calculated with a symbolic computation program, in this case, Mathematica. With MMS, we can verify that Eqs. (10) and (11) are implemented as they are written down here and that the solution converges against the prescribed analytical solution with the expected order of discretization, for our case second order in space. We have tested numerous subsystems, including just the parallel-heat-flux term isolated or just the two temperature equations isolated. In Fig. 14, the results of such a convergence test for the full system including all terms are shown starting with 8 poloidal planes and 1120 points per plane, doubling the parallel perpendicular and temporal resolution three times until we arrive at 64 poloidal planes. The supremum norm of the numerical error $Δ inf$ of all plasma quantities is converging as expected with second order. Therefore, we have strong evidence that the implementation was performed correctly. The supremum norm of the numerical error is given by $Δ inf = | f num − f MMS | ∞ / | f MMS | ∞$ with the numerical solution $f num$ obtained by GRILLIX and the analytical MMS solution $f MMS$ for each plasma field. The model was tested in simple slab geometry for the purpose of finding, e.g., typing errors in the implementation of the model. The 3D-solvers used for solving Eqs. (10) and (11) and the model with limited Braginskii closure are also tested in more complex geometries. The test performed here is not the pinnacle of verification, but it adds one additional layer of safety.