{'fraction_non_alphanumeric': 0.03831265508684863, 'fraction_numerical': 0.008476426799007444, 'mean_word_length': 4.815075608911463, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Drift reduced Landau fluid model for magnetized plasma turbulence simulations in BOUT++ framework

Recently the drift-reduced Landau fluid six-field turbulence model within the BOUT++ framework has been upgraded. In particular, this new model employs a new normalization, adds a volumetric flux-driven source option, the Landau fluid closure for parallel heat flux and a Laplacian inversion solver which is able to capture n=0 axisymmetric mode evolution in realistic tokamak configurations. These improvements substantially extended model's capability to study a wider range of tokamak edge phenomena, and are essential to build a fully self-consistent edge turbulence model capable of both transient (e.g., ELM, disruption) and transport time-scale simulations.


Introduction
Magnetized plasma is of great interest to scientists and engineers for its applications in astrophysics, space physics, material processing and fusion energy research. In this paper, we describe a fluid model that aims to study the instability, turbulence, and transport for magnetically confined plasmas. In particular, this model is designed to investigate the rich physics on the periphery of a tokamak (a toroidal device that confines high temperature plasmas with a swirling magnetic field) plasma, commonly termed as the edge plasma. The edge plasma is imperative to the success of selfsustained fusion operations as it largely controls the particle and energy exhaust, breeding and refueling, and overall fusion performance. It therefore becomes one of the most active research areas in fusion plasma and engineering nowadays.
The tokamak edge region is a complicated physics system. It normally has a radial extension of a few centimeters and consists of two distinct magnetic topologies -the closed flux region where magnetic field-lines keep circling around major axis of the torus, and the scrape-off-layer (SOL) as well as the private flux region where magnetic field-lines intersect with vessel wall and/or divertor targets. Within such a short radial distance, plasma density and temperature could drop one or two orders of magnitude, transforming from a hot, dense, nearly collisionless core plasma to a relatively cold, thin, collisional SOL plasma. The resulting strong radial inhomogeneity of plasma pressure further modifies local plasma current density profile and hence the local magnetic shear, and could excite various plasma instabilities ranging from the large size magnetohydrodynamic (MHD) instabilities to the less harmful micro-instabilities at ion and electron gyroraduis scales.
These instabilities, usually cascade into turbulence, enhance cross-field transport, alter the equilibrium profiles and degrade fusion performance. The typical edge turbulence fluctuation level may vary from O(10 −2 ) to order of unity [2] and in most cases, is largely impacted by local plasma conditions. However, global simulations suggest that fully developed edge turbulence can also be non-local, i.e., turbulence can have long-range influence on either up-gradient [3] or down-gradient [4] or both regions via turbulence spreading [5], depending on the nature of the underlying instability. In the meantime, strongly sheared flows occurred the edge region (e.g., the E × B flow in the poloidal direction and the intrinsic and/or driven rotation in the toroidal direction) are found to be crucial in terms of governing the overall edge plasma performance. Although the interaction between the sheared flows and turbulence saturation has been extensively investigated in theoretical, experimental and numerical studies, the spontaneous generation and saturation mechanisms are still not fully understood.
Last but not least, at the far-SOL and near the wall where plasma temperature is relatively cool, impurity, neutral and molecular physics become important as the atomic and molecular processes (e.g., radiation, recombination, ionization, charge-exchange, dissociation, etc) could dramatically change the local plasma dynamics and further affect upstream plasma. Therefore, tokamak edge is an inherent multi-scale, multi-physics, highly nonlinear system that requires sophisticated numerical modelling.
In recent years, a large effort has been devoted to tokamak edge modelling and significant progress has been made. Several fluid (e.g., GBS [6], TOKAM3X [7], GDB [8], GRILLIX [9]), gyro-fluid (e.g., GEMR [10] BOUT++-GLF [11]) and gyrokinetic (e.g., XGC [12], Gkeyll [13], COGENT [14]) codes have been developed to investigate different aspects of the plentiful edge related problems. However, the extremely demanding computational requirement of these high-fidelity multi-scale gyrokinetic models limits their applications in the routine global edge simulations. Concurrently, fluid models, though only contain the lowest three moments of plasma distribution functions, still capture the fundamental physics and are adequate to quantitatively describe the behavior of edge plasmas in many cases. For this reason, they are and will remain an irreplaceable component of the edge study [15].
The original Boundary Turbulence (BOUT) fluid code was first developed in the late 1990s [16,17] to explore turbulence dynamics at the tokamak boundary region (i.e., the edge). It was continuously refined over the years [18] and eventually extended to BOUT++ [1] -an open source C++ framework that provides a user friendly interface for physicists to code partial differential equations (PDEs) without worrying about numerical implementation. This successful strategy yields numerous numerical models targeted for specific research topics within the BOUT++ framework to explore a wide range of edge physics -from linear growth rate analysis to fully nonlinear, fullf turbulence simulation; from quasi-state transport to intermittent events such as ELM crashing; from particle fueling, radio frequency (RF) heating to divertor exhaust solutions. Nevertheless, the six-field turbulence model in BOUT++ framework [19] is one of the main workhorses among all BOUT++ models. Recently, it has been substantially improved -it now simulates a more complete physics model with a more practical normalization and implies a zonal field solver, a sophisticated Landau fluid parallel heat flux model and a flux-driven source option. This paper is by no means to provide a complete manual of the drift-reduced turbulence model in BOUT++, but rather a self-contained documentation that reflects current status of the model with a focus on the physics aspect to facilitate future research. The rest of the paper is organized as follows. Section 2 introduces the equations and the normalization used in the BOUT++ turbulence model. Field-aligned coordinate system that minimizes computational costs and the associated operators, as well as the boundary treatment are outlined in Section 3. Section 4 presents the latest developed Landau fluid closure on parallel heat flux for weakly collisional tokamak edge plasma in divertor configuration. The Landau fluid closure prediction is compared and contrasted with those from the widely used classical Braginskii and flux-limited models. The inversion of Laplacian operator, for both non-axisymmetric (n 0) and axisymmetric (n = 0) components, is described in Section 5. Section 6 illustrates the application of flux-driven source. Finally, Section 7 summaries the current status of the drift-reduced Landau fluid model in BOUT++ and discusses the on-going and planned efforts on future model development.

Model equations
Like the other fluid turbulence models of its kind, BOUT++ turbulence model is also based on the seminal Braginskii description of collisional, magnetized plasma [20]. In particular, BOUT++ implements the two-fluid driftreduced Braginskii equations [21,17] which eliminate high frequency dynamics (e.g., gyration and compressional Alfvén wave) by assuming that E × B and diamagnetic drifts are the lowest order flows retained in system, and therefore is suitable to investigate low frequency (ω ω ci ) turbulence and transport phenomena in the magnetized plasma. As will be discussed in Section 4 the collisionality constraint can be relaxed to some extent once the original Braginskii parallel heat flux model is replaced by the Landau fluid closure (and hence the name Landau fluid model).

Plasma equations
BOUT++'s drift-reduced Landau fluid model describes the ion scale, low frequency turbulence in a relatively collisional, magnetized, quasi-neutral (n e = Zn i ) plasma. Its equation set consists six independent primitive variables -ion density n i , electrostatic potential φ, ion parallel velocity V ,i , perturbed (parallel) magnetic flux A , electron and ion temperature T e,i , according to Here vorticity is defined as magnetic curvature total parallel current density J = J 0 + j with perturbed current density j = −∇ 2 ⊥ A /µ 0 , and ion parallel viscosity The leading ion viscosity coefficient η 0 i = 0.96nk B T i τ ii , the parallel electrical conductivity is η = F(Z)m e /(ne 2 τ e ) with F(Z) = (1 + 1.198Z + 0.222Z 2 )/(1 + 2.996Z + 0.753Z 2 ) for arbitrary Z, the E × B and ion diamagnetic drifts are Here we assume that perturbed magnetic fieldB ∇A ×b 0 is orders of magnitude smaller than the equilibrium magnetic field B 0 , a condition normally true in magnetically confined plasma, such that |B| |B 0 |, or B ≈ B 0 .
Under the drift-reduce ordering assumption, ion polarization drift is retained to leading order in order to ensure energy conservation even though its amplitude is small comparing to E × B and ion diamagnetic drifts. However, electron polarization drift along with other terms smaller by the order of m e /m i are discarded during the derivation. The electron and ion parallel heat flux q e,i will be discussed in detail in Section 4 as BOUT++ has implemented three different parallel heat flux models. The electron and ion perpendicular heat flux are often neglected in edge turbulence models since they are on the order of (ω c j τ j ) −2 1 smaller than the parallel heat flux in the collisional magnetized plasmas. Nevertheless, including cross field neoclassical thermal diffusion effects [22] in toroidal plasma simulations is available in this model. The BOUT++ Landau fluid turbulence model has a build-in option to either consistently evolve or maintain the axisymmetric plasma profiles. The formal option is called non-scale separated, or full-field model since the background plasma profiles tend to relax once the turbulence is developed and starts to continually transport particles and energy radially outward. The external particle and energy sources terms S n , S E i , S E e are hence applied when the model is used for transport time scale global full-field turbulence simulation. The latter option, on the other hand, results in the scale separated, or gradient-driven turbulence model as the gradients of the background plasma profiles which excite instabilities and turbulence are fixed over the time. These two types of model are also often refereed as fullf and δ f fluid models in a loose manner as analogous to kinetic simulations.
It is therefore reasonable and necessary to refuel particles S n and energy S E i , S E e in order to prevent the collapse of plasma profiles in transport time scale simulations. As a result, the external particle and energy flux could largely alter the local plasma profiles and further influence the plasma dynamics both in turbulence and transport time scales, so such a model is also termed as the flux-driven turbulence model.
The equation set (1)-(6) is in fact a super-set of various transport and turbulence models built upon BOUT++ framework. For instance, if one sets ∂A * /∂t = 0, it becomes a five-field electrostatic turbulence model [16,23]; if one combines density and temperature equations into total plasma pressure equation, it becomes three-field model [24]; if one makes electrostatic and axisymmetric assumptions (i.e.,∂/∂ξ = 0 where ξ is the toroidal direction) and neglected subdominant terms, it becomes 2D transport model [25].

Normalization
In the six-field model, basic normalization coefficients, e.g., the characteristic lengthL, magnetic field strengthB, densityn, electron and ion temperatureT e,i , are prescribed based on initial profiles. For most of tokamak simulations, L,B andn are the largest radius, total magnetic field strength and ion density respectively;T e,i are typically set to be 1000eV. With defined reference density and temperature, the individual pressure then can be normalized as P j =nk BT j , and the total pressureP =nk BTe ,P =P e + τP i with τ =T i /T e . Similarly,κ =Lκ,∇ =L∇, and the characteristic velocity and time areV For the electromagnetic quantities, BOUT++ chooses a normalization that minimizes the number of normalization coefficients asφ Note that in this normalization, electric field, E × B and ion diamagnetic drifts are simplŷ In several previous BOUT++ edge turbulence models [23,19], the magnetic field was normalized to background equilibrium profile B 0 instead of a fixed valueB. Consequently, electric field and E × B drift were written aŝ which require k/L B 1 (here L B is the characteristic length of background magnetic field). This condition is normally satisfied for drift-type turbulence at the tokamak edge. Nevertheless, the new normalization removes this constraint and attains simpler electric field and E × B drift expressions.
The relation between parallel current J and parallel velocities V e,i is modified aŝ To keep the parametric dependence of resistivity and ion viscosity, collision times are also normalized tot, and the resulting dimensionlessτ αβ = τ 0 αβT 3/2 α /n α where the constant coefficients (independence of n, T evolution but still radially dependent as the Coulomb logarithm λ = ln Λ is evaluated based on initial equilibrium) are The normalized (leading) ion viscosity coefficient then is and the normalized coefficient for ion heating (heat exchange) term is Similarly, the normalized Spitzer-Härm resistivity, Since the electron inertia (d 2 e ) term is small and mainly affects shear Alfvén wave, electron skin depth is assumed to be a (spatial varying) constant while the dependence of density evolution is neglected. The normalized electron skin depth square isd To further simplify equation expressions, define Here α di , α de are often associated with ion and electron diamagnetic drift related terms; while p is associated with curvature drift terms. For simplicity, the over-hat symbol for normalized quantities is omitted for the rest of the paper, and the normalized six-field equations read and In addition, to improve the numerical stability artificial forth order hyper-diffusion terms are optional on the righthand-side of the equations for each quantity. Depending on the nature of the problem and the resolution, these terms could be individually turned on to dissipate grid-scale turbulence eddies if necessary.

Coordinate, mesh, operators and boundary conditions
We briefly outline the coordinate system, mesh, operators and boundary conditions that are employed in the Landau fluid model in this section to lay the foundation for further discussion.

Coordinate and mesh
In the magnetic fusion community, there are more than a dozen choices of coordinate systems tailored for different research problems [26]; while in BOUT++, there are three coordinate systems are involved at different stages, namely the cylindrical coordinate system (R, φ, Z), the flux coordinate system (ψ, θ, ζ) and the field-aligned coordinate (x, y, z). The cylindrical coordinate system (R, φ, Z) is preferred by experimentalists due to its simplicity representing plasma equilibrium; the flux coordinate system (ψ, θ, ζ) is widely used by theorists in analytical analysis; and the field-aligned coordinate (x, y, z) is employed in the actual BOUT++ simulation due to its superior numerical efficiency. Because the cylindrical coordinate system (R, φ, Z) is mainly involved in the equilibrium calculation and mesh generation (e.g., tokamak grid generator HYPNOTOAD [27]) in BOUT++, here we focus on the latter two coordinate systems.
Before diving into the details, we need to clarify our notations. In BOUT++ notation, (R, φ, Z) is always a righthanded coordinate system as the one defined in EFIT [28] and kinetic EFIT. (ψ, θ, ζ) and (x, y, z), however, could be either right-handed or left-handed, depending on the poloidal field B θ direction. The positive toroidal ζ direction is defined as counter-clockwise when viewing from the top of torus, i.e., the opposite of the "normal" direction considered by experimentalists. The positive poloidal θ direction is defined as the direction of the polodial magnetic field B θ produced by a positive toroidal plasma current I p . In experiments, plasma current I p can be negative or positive, thus, the B θ in our notation, can also be negative or positive.
In the flux coordinate system, the tangent-basis (covariant-basis) vectors are e i = ∂r/∂i and the reciprocal-basis (contravariant-basis) vectors are e i = ∇i with i = ψ, θ, ζ. The corresponding scale factor h i = |e i | are and hence The total magnetic field in a tokamak is composed by the poloidal and toroidal magnetic fields; conveniently, it can be written as with ψ the poloidal flux and F = RB ζ . The Jacobian of flux coordinate (ψ, θ, ζ) is thus simply which again can be either positive or negative depending on the sign of poloidal field σ B θ = B θ /|B θ |, or the direction of plasma current I p .
Although the orthogonal flux coordinate (ψ, θ, ζ) works well in analytical analysis, BOUT++ employs a fieldaligned coordinate (x, y, z) to take the advantage of the strongly antiseptic feature of magnetized plasma by minimizing the resolution requirement in the field-line direction: Here ψ 0 , θ 0 are the reference poloidal flux and poloidal angle, normally set to be at the seperatrix and outboard midplane, respectively. The local pitch of the magnetic field ν is defined as With Equation (38), one can write down the transformation of the contravariant and covaraint basis (e i , e i ) between two coordinate systems, where the integrated local shear I(ψ, θ) is defined as As we therefore yield the contravariant metric tensor g i j = e i · e j = ∇i · ∇ j and the covariant metric tensor, g i j = (g i j ) −1 , Similarly, one can write down the cross products of contravariant basis of field-aligned coordinate system, Note that in the field-aligned coordinate system, constant x and z means staying on same magnetic line, as the magnetic field could simply be represented in the Clebsch form The above expression also indicates that the unit vectorb = B/|B| = B θ /(h θ |B|)e y could be in either the same or the opposite direction of e y depending on the sign of B θ . The Jacobian of (x, y, z) is given by which is equivalent to Equation (37). As mentioned previously there are two types of magnetic field lines in the tokamak edge region. Figure 1 illustrates the flux surfaces and the magnetic field lines at the tokamak edge region in a lower single null divertor configuration, the most common tokamak discharge nowadays. Inside the separatrix (white dashed), the projections of magnetic field-lines on the poloidal plane form the so-called closed flux surfaces (red lines); while outside the separatrix, the projections yield open flux surfaces as the magnetic field-lines are terminated at the divertor target plates. Depending on whether the magnetic field-lines pass through the top poloidal plane, the open flux region is further categorized into the scrape-off-layer (blue lines) and the private flux region (green lines). With these three somewhat detached regions, generating a mesh that could be parallelized efficiently is no longer straightforward. However, by carefully dividing the poloidal plane into separated blocks and connecting them based on their locations in configuration space, the simulation domain remains topologically rectangular in the index space as shown in Figure 2. To truthfully capture the geometric information, BOUT++'s tokamak meshes are normally based on ideal MHD equilibrium either attained from external solvers such as CORSICA or GATO, or reconstructed from experimental measurements via EFIT or kinetic EFIT. Auxiliary mesh generation codes (e.g., HYPERNOTOAD [27]) are developed to produce high accuracy BOUT++ meshes with desired resolution and domain. Therefore, BOUT++ is capable to simulate realistic tokamak discharges in limiter, single-null and double-null divertor configurations.
The authors remark that despite the discussion in this paper emphasizes on tokamak plasma, in principle, the generic representation of magnetic field configuration in BOUT++ makes the Landau fluid model a flexible tool to study plasmas in a variety kinds of magnetic configurations, for example, the linear devices [29], the helical plasmas and even the stellarator plasmas [30].

Operators
Here we list the most common operations involved in the fluid simulations generalized in the field-aligned coordinate system. In BOUT++, a vector v is by default in covariant form unless it is declared otherwise. However, it may be transformed in contravariant form during certain operations to reduce the overall computation.

Dot, cross, gradient, divergent, curl and Laplacian operators
The dot and cross products of two vectors u and v can be calculated either in covariant or contravariant forms, The gradient of a scalar f is most convenient to be evaluated in covariant components form (u i ) as In some cases, only the perpendicular component of ∇ f is needed which is defined as The divergence of a vector u is most convenient to be evaluated in contravariant components form (u i ) as The curl of a vector v is again most convenient to be evaluated in covariant components form (v i ) and the result is in contravariant components form (u i ), The Laplacian operator, by definition, is or, in explicit form

Poisson bracket, advection, parallel derivatives and curvature operators
In the fluid model, an important type of operator is called "Poisson bracket" defined as which physically means convection processes (e.g., E × B convection if f = φ, or magnetic flutter effects if f = A ).
Rewriting ∇ f × ∇g as and with the aid of Equations (44) and (45), the full expression for Poisson bracket becomes Because In general, an advection process of quantity f with velocity u is easier to be evaluated in contravariant form as The parallel derivative along an unperturbed magnetic fieldb = B/B is and the double derivative then becomes The derivative along the normalized perturbed magnetic fieldb =B The last type of important operators in the toroidal system is the curvature operatorb ×κ · ∇ f , whereκ is the magnetic line curvature, defined asκ =b · ∇b = −b × (∇ ×b). The curvature operator could be rewritten in a similar form of the advection operatorb where the velocity components v i =b ×κ · ∇i with i = x, y, z are calculated during mesh generation and then passed along as the prescribed inputs. Neglecting the perturbed magnetic field contribution is justified because normally it is three or more orders of magnitude smaller than the equilibrium magnetic field. It should be pointed out that in order to avoid numerical error accumulation in a sheared mesh (e.g., tokamak edge), radial derivative in BOUT++ in fact is often not done in field-aligned coordinate system but "shifted" [18] and evaluated in a so-called "quasi-ballooning" coordinate system. The implementation and advantage of shifted radial derivative is thoroughly discussed in References [1,31].

Boundary conditions
For a toroidal plasma system like tokamak, it is natural to impose periodic boundary condition in the toroidal (z) direction. In the radial direction (ψ for sheared mesh, otherwise, x), a large variety of boundary options is available within BOUT++ for different physics models, including the most commonly used Dirichlet, Neumann, Robin boundary conditions as well as the "relaxed" boundary which evolves the quantity towards the given boundary conditions at certain rate in order to prevents transient discontinuities occurring at the boundaries. In general, for the global Landau fluid model, if applicable, Dirichlet boundary is applied to the "perturbed", non-axisymmetric (i.e., n 0) quantities, and Neumann boundary is applied to the "equilibrium" axisymmetric (i.e., n = 0) quantities to minimize the fluctuations at the boundary. The parallel y direction has two distinct boundary conditions depending on its corresponding magnetic topology. As illustrated in Figure 1, inside the separatrix magnetic field-lines wind around the major axis while simultaneously swirl around minor (magnetic) axis, and nest so-called closed flux surfaces. In most cases, a magnetic field-line will not end up at the same toroidal location after finishing one poloidal turn; therefore, a twistshift operation [17,18] is used to ensure the parallel continuity in the closed flux region. In the open field-line region, i.e., the scrape-off-layer and the private flux region classical Bohm sheath criteria [32] are enforced at the boundary where magnetic field-lines are intercepted by limiter or divertor target plates, Q e,i = ±γ e,i n e,i k B v the,i T e,i , Here C s is the local sound speed, arising from the presheath electrostatic potential. In divertor configuration, theory and model predict that ion could be supersonic near the divertor plates due to the Laval nozzle effect [33,34], supersonic ion flow is also occasionally observed in UEDGE and BOUT++ transport simulations [35]. Therefore, ion parallel velocity V i is compared to the the local sound speed C s at each time step -if it is less than C s , then |V i | = C s ; otherwise, it is free evolving, i.e., supersonic flow is allowed. The electron/ion sheath energy transmission factors γ e,i depend on secondary emission coefficient, local electron and ion temperature ratio, target plate floating potential [36].
The typical values used in BOUT++ simulations are γ e 7, γ i 2.5. As will be discussed in next section, the implementation of this boundary condition depends on the specific parallel heat flux model to be used. In addition the following boundary conditions on density, vorticity and electrostatic potential are applied to eliminate the potential build-up of large discontinuities at the boundary,

Braginskii (Spitzer-Härm) heat flux model
The Braginskii transport model [20] derived under the assumption that plasma is magnetized and collsional requires that L ⊥ ρ and L λ mfp . Here ρ is the particle (electron and ion) gyroradius, λ mfp is the particle mean free path, and L ⊥ , L are the characteristic lengths in perpendicular and parallel directions. Under the strong collisional limit, the well-known Braginskii (or, Spitzer-Härm) heat flux model is deduced by assuming the perturbed particle distribution slightly deviates from the equilibrium Maxwellian distribution, which implies that heat flux depends on local temperature gradient, and heat (or, energy) is transported solely through collisions (or, random walk process) among particles.
Recent particle simulations indicate that Braginskii parallel heat flux model is reasonably accurate (i.e., less than 10% error) when the parallel inhomogeneity length is at least on the two order of magnitude longer than the particle mean free path (i.e., L /λ mfp ≥ O(10 2 )) and becomes invalid when L /λ mfp ≤ O(10) [37].
In the tokamak edge region, L ⊥ is the equilibrium length such as L n , L T which are typically a few centimeters, while L ∼ 2πqR in large aspect ratio limit is on the order of tens meters where q ≈ 4 ∼ 6 is the safety factor and R is the major radius of the machine. Assuming that plasma temperature is isotropic (i.e., T ⊥ = T = T ), then ρ = v th /ω c ∝ m 1/2 T 1/2 and λ mfp ∼ v th τ ∝ T 2 /n is independent of mass. As modern and future tokamaks produce higher and higher temperature plasmas, the latter requirement (L λ mfp ) is no long be satisfied due to the strong temperature dependence of λ mfp . For example, Table 1

Flux-limited model
Intuitively, in a collisionless plasma, parallel heat transport is mainly done by particles freely moving along the magnetic field-lines which also sets a upper limit of heat flux. In fluid description, this free streaming heat flux is defined as q FS = nk B T v th where v th is the thermal speed. Hence, a flux-limited expression for parallel heat flux is constructed to constrain the unrealistically large value predicted by Braginskii model [40], where α ∈ [0.1, 1] is a parameter to be adjusted in different scenarios. However, Equation (70) becomes singular at where q B + αq FS ≈ 0 as by definition, q FS is always positive while q B can be negative depending on local temperature gradient and |q B | q FS in the weakly collisional regime. In practice, q FS is replaced with the "diffusive" free streaming heat flux q FS(d) = −nv th L ∇ (k B T ) which ensures it has the same sign as q B at all locations. Consequently, the resulting flux-limited parallel heat flux is again dependent on local temperature gradient, similar to the Braginskii model. It can be rewritten as q FL j = −κ j eff ∇ (k B T j ) with the effective thermal conductivity

Landau fluid closure
Although the flux-limited parallel heat flux model overcomes the overestimation of Braginskii model, it is still an ad hoc model. A more rigorous approach is to develop a theory-based parallel heat flux model, for example, the Landau fluid model for weakly collisional plasmas [41] where α e 0.499 and α i 0.409. This Landau fluid closure is similar to the closure by incorporating the Krook collision operator [42,43] with slightly different coefficients and an extension of Hammett-Perkins closure [44] -the first and the simplest Landau fluid closure that describes plasma kinetic responses in the collsionless limit, As a result, it seamlessly bridges collsional Braginskii expression (69) and collsionless Hammett-Perkins closure (73) as shown in Figure 3. Comparing with the exact kinetic closures with collisional effects [45], this closure predicts a slightly lower heat flux (with a maximum difference by 25%) in the intermediate collisionality regime. Implementing such a closure in BOUT++ turbulence model is numerically challenging. This is because Landau fluid closures are usually entail wave-vector space (i.e., k−space) and implicitly assume periodicity along magnetic field-lines; while BOUT++ turbulence model calculates quantities in configuration space. A novel solution to resolve these issues is the fast non-Fourier method [46] which approximates parallel heat flux with summation of Lorentizians as where Here ξ j = λ m f p, j /α j and (α n , β n ) are optimized fitting coefficients of as the ones given in Table 2 of Reference [47]. It is therefore straightforward to transform the above Lorentizians back to configuration space by replacing ik and k 2 with ∇ and −∇ 2 , and the evaluation of parallel heat flux in k−space becomes solving a series of boundary value problems in configuration space along unperturbed magnetic field-lines. The fast non-Fourier method has been developed and validated with a one-dimensional uniform grid [47]. As illustrated in Figure 4, it attains more accurate approximation and reaches further into the collisionless regime by adding more Lorentizians. In most tokamak edge simulations, N = 7 is sufficient as it is valid for ξk ∼ k λ m f p < 10 4 with a relative error about 2%. One can further push the applicable limit to k λ m f p < 10 7 and/or reduce the relative error to about 1% by using N = 12 Lorentizians. This treatment of Landau fluid closure has now been extended and implemented in BOUT++ transport and turbulence codes to deal with three-dimensional divertor geometry. To highlight the discrepancies between the original Braginskii heat flux model, flux-limited model and the Landau fluid closure in a realistic tokamak geometry, two scenarios are considered -a relatively collisional C-Mod discharge ( Figure 5) and a less collisional DIII-D discharge ( Figure 6) based on their typical H-mode parameters. In both scenarios, electron parallel heat flux is evaluated assuming there is a relatively large (30%) temperature perturbation is on the top portion of the equilibrium profiles (e.g., to imitate vertical displacement event-type instability).
The first key observation is that Braginskii model predicts a parallel heat flux only at the tokamak top region where temperature gradient is finite (thus, "local" model) and its amplitude is about two orders of magnitude larger than either flux-limited or the Landau fluid estimate. The flux-limited model seems able to effectively scale down the heat flux amplitude close to that predicted by the Landau fluid closure (Figures 5(e) versus (f), Figures 6(e) versus (f)); nevertheless, it is still local. On the other hand, the Landau fluid closure obviously exhibits the "nonlocal" effect as the predicted parallel heat flux extends to the region where local temperature gradient vanishes (e.g., the tokamak bottom region). More interestingly, our tests also demonstrate the transition from local-dominant transport to nonlocal-dominant transport. In C-Mod plasma, parallel heat flux still marginally depends on the local temperature gradient (Figures 5(f)). As the plasma becomes less collisional, this local effect weakens and the nonlocal effect dominates as expected (Figures 5(f) versus 6(f)).
Among the three parallel heat flux models, only the Landau fluid closure is able to consistently enforce the ambipolar sheath heat flux boundary condition q sheath at the divertor targets (i.e., Equation 66) as shown in Figures 5(i) and 6(i). This is again due to the distinct inherent transport characteristics. Because Braginskii and flux-limited models are local, q sheath is actually not enforced on parallel heat flux but on the plasma temperature as ∇ T = ±q sheath /κ which tends to produce a steep parallel gradient at the divertor targets. While for Landau fluid closure, q sheath is provided directly as the Dirichlet boundary condition and hence results in a smoother parallel heat flux profile.
Given the dramatic discrepancies between three parallel heat flux models, one would expect that plasma parallel dynamics may also be altered significantly. Perhaps a more important quantity is the parallel gradient of parallel heat flux (∇ q ) since it directly influences temperature evolution. Our initial analysis suggests that Landau fluid closure results in a nearly one order of magnitude weaker heat conduction comparing to the value given by the flux-limited model both in the closed flux region and in the SOL, as attested by the substantial diminishing amplitudes of ∇ q (red dotted lines) in Figure 7. However, the impact of the improved Landau fluid closure on edge instability, turbulence, transport and most importantly, divertor heat flux width [48] doesn't fit in the scope of this paper and will be left for a future study.
Although the Landau fluid closure is arguably the most sophisticated treatment of parallel heat flux in nowadays fluid edge turbulence models, it is worth to clarify that there are still a few subtleties in this approach. For instance, the temperature anisotropy effect is neglected in order to be consistent with Braginskii transport equations, and the classical Bohm's sheath criteria for collisional plasma is used at the boundaries in the open field-line region. The former can be addressed in a more complete BOUT++ gyro-Landau-fluid turbulence model. Although one may argue that this Bohm's sheath boundary condition is only valid in collisional fluid regime, the mean kinetic (from the half-Maxwellian) and fluid velocities are actually very similar (v = √ 8T/m i π vs c s ) because ions that reach the divertor plate will escape. While the heat flux predicted by Landau fluid closure (72) is found to have a good agreement with the result from a Fokker-Planck code with a small amplitude perturbation in weakly collisional regime [41], a further validation with extended parameter regime will be performed and reported in future publications.

Normalization
For the sake of completeness, the normalization of different parallel heat flux models are listed below. In BOUT++ transport and turbulence codes, q j is naturally normalized withq j =nk BVAT j and the normalized parallel heat flux for Braginskii expression and flux-limited model along with their thermal conductivities are with κ e = 3.2τ 0 ei T 3/2 e v 2 th,e , κ e = 3.9τ 0 ii T 3/2 i v 2 th,i , κ j eff = αn j v th, j L κ j αn j v th, j L + κ j . Figure 7: q and ∇ q predicted by three models along field-lines (a) at pedestal top (ψ N = 0.95) and (b) in the SOL (ψ N = 1.01) in the DIII-D test (i.e., Figure 6).
Here the plasma density and temperature n j , T j , thermal speed v th, j and parallel characteristic length L are normalized ton,T j ,V A andL respectively. For the Landau fluid closure, Equation (77) that solves Lorentizians of parallel heat flux after normalization becomes where n j , T j , ξ j and v th, j are again normalized ton,T j ,L andV A .

Field solver
In the drift-reduced fluid formalism, Laplacian inversion is required to evaluate electrostatic potential φ from general vorticity , and parallel magnetic potential A (or perturbed magnetic flux ψ=−A ) from modified magnetic flux A * if electron inertia is included in the model, Therefore, Laplacian inversion solvers are designed in BOUT++ to solve both equations in a generalized form Since BOUT++ employs the field aligned coordinate system, the perpendicular component of Laplacian operator used in the vorticity equation and Ohm's law, by definition, is where G i = ∇ 2 i for i = x, y, z and the ∂ 2 xy term is omitted as g xy = 0. Instead of directly solving in the matrix inversion problem within 3D configuration space, BOUT++ separates the unknown quantify f into non-axisymmetric (n 0) and axisymmetric (n = 0) components to improve the numerical performance by reducing the dimensionality of the problem. Here n represents the toroidal mode number.

Non-axisymmetric (n 0) field
For the non-axisymmetric field, considering the strongly anisotropic nature of magnetized plasmas, the parallel derivatives are normally several orders of magnitude smaller than the perpendicular derivatives, i.e., k k ⊥ . Therefore we may neglect any ∂/∂y terms and also the G x , G z terms (as the grid is normally uniform in x and z direction), and Equation (84) becomes (85) Rewriting above equation in flux coordinate system to drop all the I terms which may introduce large error in a sheared mesh, and then applying Fourier transform along the toroidal direction, yields a 1D boundary value problem and can be solved efficiently with the tridiagonal solvers [31].

Axisymmetric (n = 0) field
However, the above approach doesn't apply to the axisymmetric k z = 0 (or n = 0) zonal field component. Intuitively, this means that the compressed 1D profile (x distribution) does not contain enough information to reconstruct a 2D structure (xy−plane); therefore, one needs to include y variation in order to correctly solving the zonal field component. In other similar edge turbulence codes, this process is normally done with a hard-coded algebraic multigrid solver [49,50,51]. BOUT++ utilizes external libraries (e.g., PETSc [52], HYPRE [53]) to solve this problem thus has a more flexible option on the numerical algorithm to be used [54,31]. In BOUT++ 2D transport code, the relaxation method is used to solve the zonal potential equation for a steady-state solution. [55] For n = 0 mode, ∂/∂z terms vanishes, and Equation (84) becomes Above equation is then discretized either based on the second order finite difference scheme (thus, 5-point stencil) or the forth order finite difference scheme (thus, 9-point stencil) and solved as a 2D boundary value problem.
To verify that the solver is implemented properly, an inverse-and-forward test on electrostatic potential is performed. That is, giving an self-consistent set of profiles (n i , T i , B) and assuming = 0, evaluate the right-hand-side of Equation (81), i.e., and calculate φ by solving the 2D boundary value problem; then re-evaluate the left-hand-side of Equation (81) (denoted as RHS ) based on φ; and finally the absolute error δ RHS = RHS − RHS is used to quantify the performance of the solver in this quasi-realistic setting. In this test, Equation (81) is discretized with a 5-point stencil, and the generalized minimal residual (GMRES) method and the successive over-relaxation (SOR) method in the PETSc library are picked as the matrix solver and the preconditioner. Dirichlet boundary condition (φ = constant) is applied on the outside boundary (wall-side) and Neumann boundary condition (∂φ/∂ψ = 0 ∝ V E,p ) is enforced on the inside boundary (core-side) to eliminate any external poloidal angular momentum sourcing. Additional options, such as the absolute and relative error tolerances A , R set the iteration criteria, could be used to adjust the overall solver accuracy. Figure 8 shows a typical test in a shift-circular configuration with resolution (n x , n y ) = (516, 64). The absolute error δ RHS ∼ 10 −2 (Figure 8(f)) is about eight orders of magnitude smaller than the initial input (Figure 8(c)), suggesting that the axisymmetric field solver is correctly implemented. Moreover, the resulting error δ RHS tends to collocate at where RHS is finite. A convergence test is carried out by varying the relative error tolerance R with and without the ∇ ⊥ n i · ∇ ⊥ φ cross term on the left-hand-side of Equation (81), and the result is summarized in Figure 10(a). Not surprisingly, the axisymmetric field solver converges. As R decreases, the maximum absolute error of the test |δ RHS | max decreases linearly as well till it saturates at 10 −2 level. As expected, the solver slows down as R becomes smaller, i.e., the number of iterations required to achieve certain criterion increases. However, the apparent discrepancy between the two red lines indicates that the number of iterations is reduced substantially when the cross term is neglected (i.e., applying Boussinessq approximation -a widely used but not fully justified assumption in reduced-fluid research) especially with a fixed relatively small R . It is worth to point out that the exact number of iterations depends on quite a few factors, for example, grid resolution, smoothness and consistency of the profiles (RHS and n i ), choices of matrix solver and preconditioner and so on.
Same inverse-and-forward test is also performed in a C-Mod-like lower single null configuration with resolution (n x , n y ) = (68, 128) as illustrated in Figures 9 and 10(b). Similar to the shift-circular case, the axisymmetric field solver is able to calculate electrostatic potential φ and reproduce the right-hand-side of Equation (81) in a good precision although in this case, error seems to concentrate at the far private flux region (likely due to the choices of initial profiles and boundary conditions). Interestingly, in the convergence test, no significant discrepancy of the number of iterations with and without the cross term is observed. The new axisymmetric field solver is further benchmarked against relaxation method [55] currently implemented in BOUT++'s transport model. As plotted in Figure 11, for a given initial potential and vorticity profiles φ 0 , 0 , both methods are able to recover the initial potential profile. The overall small discrepancy between two solutions is likely due to the slightly shifted location of boundary condition implementation -while relaxation method enforces boundary condition at the boundary grid points, the axisymmetric field solver enforces boundary condition at the boundary cell centers.
Authors would like to point out that the new field solver is able to invert the a more accurate form of vorticity, e.g., including the n = 0 component of the cross term contributed from axisymmetric ion density and electrostatic potential in Equation (29). For years, edge turbulence models either drops the crossing term ∇ ⊥ φ · ∇ ⊥ n i and/or the ion density dependence in front of ∇ 2 ⊥ φ term (i.e., so called Boussinesq approximation) in order to simplify the expression of vorticity, hence the numerical calculation, However, the validity of Boussinesq approximation in global edge turbulence simulations is never rigorously justified,  especially considering that radial density profile n may have a steep gradient at the pedestal region. Recent studies suggest that Boussinesq approximation has limited influence (∼ 10%) on the overall turbulent fluctuation level with [56] and without [57] hot ion in a circular geometry. Our preliminary tests indicate that the cross term and the 1/n i dependence do influence the inverted axisymmetric electrostatic potential and the corresponding radial electric field (i.e., zonal field) as shown in Figure 12. A through study of the impact of Boussinesq approximation on edge instability, turbulence, flow generation and saturation is subject of future research.

Flux driven source
In the edge turbulence simulations, plasma is constantly transported outward once instability is excited, it is therefore reasonable to refuel particles and energy from the core side in order to prevent the collapse of density and temperature profiles in the transport time scale simulations. Meanwhile, plasma sources are also necessary to simulate certain transient events, e.g., the thermal quench phase of tokamak disruption when the hot core plasma is pushed towards edge in a very short time-scale ( 100 µs).
The newly added sourcing options therefore largely expands the range of physics problems that BOUT++ edge model is able to investigate. A good demonstration of this new capability is the preliminary study of divertor heat load response to excessive heating which is an important issue in tokamak disruption mitigation. A nonlinear electromagnetic simulation is first carried out in a lower single null DIII-D H-mode plasma without sourcing for the first 0.14 ms. For the next 80 µs, 1 GW power, or 80 kJ energy in total (roughly 15% of the total stored plasma thermal energy) is injected to the pedestal top region (0.85 < ψ N < 0.96) to mimic intensive energy flux coming from the core in the thermal quench phase. For simplicity, heating power is equally partitioned between electrons and ions, and assumed to deposit uniformly in the poloidal and toroidal directions while take a Gaussian shape in the radial direction. At t = 0.22 ms, energy source is turned off. Flux-limited parallel heat flux model is used in this particular simulation. Figure 13 plots the time history of outboard midplane electron temperature T e,OMP at the pedestal top (ψ N = 0.95) and the corresponding peak magnitude of the heat flux at outer divertor target plate. A monotonic, almost linear increasing of T e,OMP is observed once the power injection starts, followed by a much slower decreasing phase once the source is turned off. The divertor heat flux, on the other hand, takes a while to response to the change of pedestal temperature. The 30 µs lag (e.g., between the peak temperature and peak divertor heat load) is mainly due to the parallel conduction process where the plasma inside the separatrix first enters the SOL primarily near the outboard midplane and then is transported to the divertor target. In this run, doubled pedestal top temperature yields more than 40 times larger heat flux and doubled heat flux width. The dramatic increasing of heat load, as well as the broadening of heat flux width, are caused by the enhanced turbulence activity at the tokamak edge as illustrated in Figure 14. Before the intensive energy injection, the system is marginally stable to the intermediate-n peelingballooning modes localized at the steep gradient region right inside the separatrix (Figure 14(a)). Once the heating is on, elevated temperature profiles steepen edge plasma pressure gradients and destabilize these previously marginal stable intermediate-n modes. Despite these modes are strongly sheared by the edge E r field, they are still able to spread both inward to the pedestal top and outward across the separatrix into the scrape-off-layer, and therefore, amplifies cross-field heat transport (Figure 14(b)). As the pedestal plasma is substantially heated, low-n modes are eventually excited near the pedestal top, providing a more violent yet effective turbulent transport channel and resulting an even larger heat load on the divertor target plate (Figure 14(c)). Although the details of the underlying physics in process is yet to be explored, the versatile flux-driven source option provides a unique opportunity to study edge plasma dynamics in this pulsed over-driven setting.

Summary
The drift-reduced two fluid six-field turbulence model in BOUT++ framework is perhaps the most widely used tokamak edge turbulence code in the magnetic fusion community. It evolves plasma density, electron and ion temperature, ion parallel velocity, electrostatic potential and parallel magnetic flux in a realistic tokamak magnetic configuration. In this paper, we present in detail the current status of this model in a self-contained manner, including the model equations, normalization, coordinate system, operators, boundary conditions and the recently added features such as flux-driven source, Landau fluid closure on parallel heat flux and the axisymmetric field solver. These new features greatly extend the model capability to simulate transport time scale plasma dynamics (i.e., shear flow generation, turbulence saturation and profile evolution) and to study critical edge problems such as divertor heat flux solutions.
It should be noted that BOUT++ remains an active project. Numerical and physics models within BOUT++ are continuously refined. For example, in the numerical side, GPU-capability and a new 3D field solver are under development; while on the physics side, a comprehensive fluid neutral model and a ML-enabled surrogate model [58] of an even more sophisticated kinetic Landau fluid closure [59] are also on the way.