Drift reduced Landau fluid model for magnetized plasma turbulence simulations in BOUT++ framework
DDrift reduced Landau fluid model for magnetized plasma turbulence simulationsin BOUT ++ framework Ben Zhu a, ∗ , Haruki Seto b , Xue-qiao Xu a , Masatoshi Yagi b a Lawrence Livermore National Laboratory, Livermore, California 94550, USA b National Institutes for Quantum and Radiological Science and Technology, Rokkasho, Aomori 039-3212, Japan
Abstract
Recently the drift-reduced Landau fluid six-field turbulence model within the BOUT ++ framework [1] has beenupgraded. In particular, this new model employs a new normalization, adds a volumetric flux-driven source option, theLandau fluid closure for parallel heat flux and a Laplacian inversion solver which is able to capture n = Keywords: tokamak edge; turbulence; transport; Landau fluid; Braginskii equations.
1. 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 instabil-ity, turbulence, and transport for magnetically confined plasmas. In particular, this model is designed to investigate therich physics on the periphery of a tokamak (a toroidal device that confines high temperature plasmas with a swirlingmagnetic field) plasma, commonly termed as the edge plasma. The edge plasma is imperative to the success of self-sustained fusion operations as it largely controls the particle and energy exhaust, breeding and refueling, and overallfusion performance. It therefore becomes one of the most active research areas in fusion plasma and engineeringnowadays.The tokamak edge region is a complicated physics system. It normally has a radial extension of a few centimetersand consists of two distinct magnetic topologies – the closed flux region where magnetic field-lines keep circlingaround major axis of the torus, and the scrape-o ff -layer (SOL) as well as the private flux region where magneticfield-lines intersect with vessel wall and / or divertor targets. Within such a short radial distance, plasma density andtemperature could drop one or two orders of magnitude, transforming from a hot, dense, nearly collisionless coreplasma to a relatively cold, thin, collisional SOL plasma. The resulting strong radial inhomogeneity of plasma pressurefurther modifies local plasma current density profile and hence the local magnetic shear, and could excite variousplasma 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 fluctuationlevel may varies from O (10 − ) to order of unity 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 canhas long-range influence on either up-gradient or down-gradient or both regions via turbulence spreading, dependingon the nature of the underlying instability. In the meantime, strongly sheared flows at 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 ∗ Corresponding author.
E-mail address: [email protected]
Preprint submitted to February 10, 2021 a r X i v : . [ phy s i c s . p l a s m - ph ] F e b n terms of instability suppression and turbulence saturation, although the spontaneous generation of these shearedflows, likely linked back to the turbulence, is still not fully understood. Moreover, magnetic geometry or plasmashaping (e.g., elongation, triangularity) could also impact the edge instability. Last but not least, at the far-SOL andnear the wall where plasma temperature is relatively cool, impurity, neutral and molecular physics become importantas the atomic and molecular processes (e.g., radiation, recombination, ionization, charge-exchange, dissociation, etc)could dramatically change the local plasma dynamics and further a ff ect upstream plasma. Therefore, tokamak edge isan inherent multi-scale, multi-physics, highly nonlinear system that requires sophisticated numerical modelling.In recent years, a large e ff ort has been devoted to tokamak edge modelling and significant progress has beenmade. Several fluid (e.g., GBS [2], TOKAM3X [3], GDB [4], GRILLIX [5]), gyro-fluid (e.g., GEMR [6] BOUT ++ -GLF [7]) and gyrokinetic (e.g., XGC [8], Gkeyll [9],COGENT [10]) codes have been developed to investigate di ff erentaspects of the plentiful edge related problems. However, the extremely demanding computational requirement ofthese 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, stillcapture the fundamental physics and are adequate to quantitatively describe the behavior of edge plasmas in manycases. For this reason, they are and will remain an irreplaceable component of the edge study.The original Boundary Turbulence (BOUT) fluid code was first developed in the late 1990s [11, 12] to exploreturbulence dynamics at the tokamak boundary region (i.e., the edge). It was continuously refined over the years [13]and eventually extended to BOUT ++ [1] – an open source C ++ framework that provides a user friendly interfacefor physicists to code partial di ff erential equations (PDEs) without worrying about numerical implementation. Thissuccessful strategy yields numerous numerical models targeted for specific research topics within the BOUT ++ frame-work to explore a wide range of edge physics – from linear growth rate analysis to fully nonlinear, full- f turbulencesimulation; from quasi-state transport to dynamic ELM crashing; from particle fueling, radio frequency (RF) heatingto particle and energy loss at wall and divertor plates. Nevertheless, the six-field turbulence model in BOUT ++ frame-work [14] 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 accurate 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 ++ , butrather a self-contained documentation that reflects current status of the model with a focus on the physics aspect tofacilitate future research. The rest of the paper is organized as follows. Section 2 introduces the equations and thenormalization used in the BOUT ++ turbulence model. Field-aligned coordinate system that minimizes computationalcosts and the associated operators, as well as the boundary treatment are outlined in Section 3. Section 4 presentsthe latest developed Landau fluid closure on parallel heat flux for weakly collisional tokamak edge plasma in divertorconfiguration. The Landau fluid closure prediction is compared and contrasted with those from the widely usedclassical Braginskii and flux-limited models. The inversion of Laplacian operator, for both non-axisymmetric ( n (cid:44) n =
0) components, is described in Section 5. Section 6 illustrates the application of flux-drivensource. Finally, Section 7 summaries the current status of the drift-reduced Landau fluid model in BOUT ++ anddiscusses the on-going and planned e ff orts on future model development.
2. Model equations
Like the other fluid turbulence models of its kind, BOUT ++ turbulence model is also based on the seminal Bra-ginskii description of collisional, magnetized plasma [15]. In particular, BOUT ++ implements the two-fluid drift-reduced Braginskii equations [16, 12] which eliminate high frequency dynamics (e.g., gyration and compressionalAlfv´en wave) by assuming that E × B and diamagnetic drifts are the lowest order flows retained in system, andtherefore is suitable to investigate low frequency ( ω (cid:28) ω ci ) turbulence and transport phenomena in the magnetizedplasma. As will be discussed in Section 4 the collisionality constraint can be relaxed to some extent once the originalBraginskii parallel heat flux model is replaced by the Landau fluid closure (and hence the name Landau fluid model). BOUT ++ ’s drift-reduced Landau fluid model describes the ion scale, low frequency turbulence in a relativelycollisional, magnetized, quasi-neutral ( n e = Zn i ) plasma. Its equation set consists six independent primitive variables2 ion density n i , electrostatic potential φ , ion parallel velocity V (cid:107) , i , perturbed (parallel) magnetic flux A (cid:107) , electron andion temperature T e , i , according to ∂∂ t n i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ n i − n i B ˆ b × κ · ∇ ⊥ φ − ZeB ˆ b × κ · ∇ ⊥ P i − n i B ∇ (cid:107) (cid:32) V (cid:107) i B (cid:33) + S n , (1) ∂∂ t (cid:36) = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ (cid:36) + B ∇ (cid:107) (cid:32) J (cid:107) B (cid:33) + ˆ b × κ · ∇ P − ˆ b × κ · ∇ π ci − Ω i (cid:104) n i Ze V D i · ∇ (cid:16) ∇ ⊥ φ (cid:17) − m i Ω i ˆ b × ∇ n i · ∇ V E (cid:105) + Ω i (cid:104) V E · ∇ (cid:16) ∇ ⊥ P i (cid:17) − ∇ ⊥ ( V E · ∇ P i ) (cid:105) , (2) ∂∂ t V (cid:107) i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ V (cid:107) i − m i n i ∇ (cid:107) P − m i n i B ∇ (cid:107) (cid:32) π ci B (cid:33) − V D i · ∇ V (cid:107) i − V (cid:107) i S n n i , (3) ∂∂ t A ∗(cid:107) = −∇ (cid:107) φ + η (cid:107) µ ∇ ⊥ A (cid:107) + en e ∇ (cid:107) P e + . k B e ∇ (cid:107) T e , (4) ∂∂ t T i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ T i + n i k B ∇ (cid:107) q (cid:107) i + m e m i Z τ e ( T e − T i ) − T i (cid:34)(cid:32) B ˆ b × κ (cid:33) · (cid:32) ∇ φ + Zen i ∇ P i + k B Ze ∇ T i (cid:33) + B ∇ (cid:107) (cid:32) V (cid:107) i B (cid:33)(cid:35) − π ci k B n i (cid:34) √ B ∇ (cid:107) (cid:16) √ BV (cid:107) i (cid:17) − k B Zen i B ˆ b · ∇ n i × ∇ T i (cid:35) − Ω i T i V (cid:107) i ˆ b × κ · ∇ V (cid:107) i + S Ei n i − T i S n n i , (5) ∂∂ t T e = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) e ˆ b (cid:33) · ∇ T e + n e k B ∇ (cid:107) q (cid:107) e + .
71 2 T e en e B ∇ (cid:107) (cid:32) J (cid:107) B (cid:33) − m e m i τ e ( T e − T i ) − T e (cid:34)(cid:32) B ˆ b × κ (cid:33) · (cid:32) ∇ φ − en e ∇ P e − k B e ∇ T e (cid:33) + B ∇ (cid:107) (cid:32) V (cid:107) e B (cid:33)(cid:35) + n e k B η (cid:107) J (cid:107) + S Ee n e − T e S n n e . (6)Here vorticity (cid:36) is defined as (cid:36) = n i m i B (cid:32) ∇ ⊥ φ + n i ∇ ⊥ φ · ∇ ⊥ n i + Zen i ∇ ⊥ P i (cid:33) , (7)magnetic curvature κ = ˆ b · ∇ ˆ b = ∇ × ˆ b × ˆ b with ˆ b = B / B , (8)modified magnetic flux A ∗(cid:107) = (cid:16) − d e ∇ ⊥ (cid:17) A (cid:107) with d e = c /ω pe , (9)total parallel current density J (cid:107) = J (cid:107) + j (cid:107) with perturbed current density j (cid:107) = −∇ ⊥ A (cid:107) /µ , and ion parallel viscosity π ci = η i (cid:34) κ · ( V E + V D i ) − √ B ∂ (cid:107) ( √ BV (cid:107) i ) (cid:35) . (10)The leading ion viscosity coe ffi cient η i = . nk B T i τ ii , the parallel electrical conductivity is η (cid:107) = F ( Z ) m e / ( ne τ e )with F ( Z ) = (1 + . Z + . Z ) / (1 + . Z + . Z ) for arbitrary Z , the E × B and ion diamagnetic drifts are V E = ˆ b × ∇ ⊥ φ/ B , V D i = ˆ b × ∇ ⊥ P i / ( Zen i B ) . (11)Under the drift-reduce ordering assumption, the ion polarization drift is also retained to leading order in orderto ensure energy conservation even though its amplitude is small comparing to E × B and ion diamagnetic drifts.The electron polarization drift along with other terms smaller by the order of m e / m i , however, are discarded duringthe derivation. The electron and ion parallel heat flux q (cid:107) e , i will be discussed in detail in Section 4 as BOUT ++ has implemented three di ff erent parallel heat flux models. The electron and ion perpendicular heat flux are oftenneglected in edge turbulence models since they are on the order of ( ω c j τ j ) − (cid:28) ff usion e ff ects [17]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 callednon-scale separated, or full-field model since the background plasma profiles tend to relax once the turbulence isdeveloped and starts to continually transport particles and energy radially outward. The external particle and energysources terms S n , S Ei , S Ee are hence applied when the model is used for transport time scale global full-field turbulencesimulation. The latter option, on the other hand, results in the scale separated, or gradient-driven turbulence modelas 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 full- f and δ f fluid models in a loose manner as analogous tokinetic simulations.It is therefore reasonable and necessary to refuel particles S n and energy S Ei , S Ee in order to prevent the collapseof plasma profiles in transport time scale simulations. As a result, the external particle and energy flux could largelyalter 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 ∗(cid:107) /∂ t =
0, it becomes a five-field electrostatic turbulence model [11, 18]; if onecombines density and temperature equations into total plasma pressure equation, it becomes three-field model [19]; ifone makes electrostatic and axisymmetric assumptions (i.e., ∂/∂ξ = ξ is the toroidal direction) and neglectedsubdominant terms, it becomes 2D transport model [20]. In the six-field model, basic normalization coe ffi cients, e.g., the characteristic length ¯ L , magnetic field strength ¯ B ,density ¯ n , electron and ion temperature ¯ T e , i , are prescribed based on initial profiles. For most of tokamak simulations,¯ L , ¯ B and ¯ n are the largest radius, total magnetic field strength and ion density respectively; ¯ T e , i are typically setto be 1000 eV . With defined reference density and temperature, the individual pressure then can be normalized as¯ P j = ¯ nk B ¯ T j , and the total pressure ¯ P = ¯ nk B ¯ T e , ˆ P = ˆ P e + τ ˆ P i with τ = ¯ T i / ¯ T e . Similarly, ˆ κ = ¯ L κ , ˆ ∇ = ¯ L ∇ , and thecharacteristic velocity and time are ¯ V A = ¯ B / (cid:112) µ m i ¯ n , ¯ t = ¯ L / ¯ V A . (12)For the electromagnetic quantities, BOUT ++ chooses a normalization that minimizes the number of normalizationcoe ffi cients as ¯ φ = ¯ L ¯ B / ¯ t , ¯ E = ¯ B ¯ V A , ¯ (cid:36) = m i ¯ n / ¯ t , ¯ J = ¯ B / ( µ ¯ L ) , ¯ A (cid:107) = ¯ L ¯ B , ¯ η = µ ¯ L ¯ V A . (13)Note that in this normalization, electric field, E × B and ion diamagnetic drifts are simply ˆ E = − ˆ ∇ ˆ φ, ˆ V E = ˆ b × ˆ ∇ ˆ φ ˆ B , ˆ V D i = k B ¯ T i Ze ¯ B ¯ L ¯ V A ˆ b × ˆ ∇ ˆ P i ˆ n i ˆ B . (14)In the previous BOUT ++ edge turbulence model, the magnetic field was normalized to background equilibrium profile B instead of a fixed value ¯ B . Consequently, electric field and E × B drift were written as ˆ E = − ˆ ∇ (cid:16) B ˆ φ (cid:17) B (cid:39) − ˆ ∇ ˆ φ, ˆ V E = ˆ b × ˆ ∇ (cid:16) B ˆ φ (cid:17) B (cid:39) ˆ b × ˆ ∇ ˆ φ (15)which require k / L B (cid:28) L B is the characteristic length of background magnetic field). This condition is normallysatisfied for drift-type turbulence at the tokamak edge. Nevertheless, the new normalization removes this constrainand attains more accurate electric field and E × B drift.The relation between parallel current J (cid:107) and parallel velocities V (cid:107) e , i is modified asˆ V (cid:107) e = ˆ V (cid:107) i − ¯ B µ ¯ L ¯ n ¯ V A e ˆ J (cid:107) ˆ n e . (16)4o keep the parametric dependence of resistivity and ion viscosity, collision times are also normalized to ¯ t , andthe resulting dimensionless ˆ τ αβ = τ αβ ˆ T / α / ˆ n α where the constant coe ffi cients (independence of n , T evolution but stillradially dependent as the Coulomb logarithm ln Λ is evaluated based on initial equilibrium) are τ ee = . × ¯ T / e λ ¯ n ¯ t , τ ei = . × ¯ T / e Z λ ¯ n ¯ t , τ ii = . × √ µ ¯ T / i Z λ ¯ n ¯ t . (17)The normalized (leading) ion viscosity coe ffi cient then is η i = . n i ˆ T i τ ˆ τ ii = (cid:15) g ˆ T / i where (cid:15) g = . τ ii τ i , (18)and the normalized coe ffi cient for ion heating (heat exchange) term is2 m e m i ˆ τ ei = (cid:15) q ˆ n e ˆ T / e where (cid:15) q = m e m i τ ei . (19)Similarly, the normalized Spitzer-H¨arm resistivity, η (cid:107) = η (cid:107) ˆ T / e where η (cid:107) = . × − F ( Z ) Z λ ¯ T / e µ ¯ L ¯ V A . (20)Since the electron inertia ( d e ) term is small and mainly a ff ects shear Alfv´en wave, electron skin depth is assumed tobe a (spatial varying) constant while the dependence of density evolution is neglected. The normalized electron skindepth square is ˆ d e = c ω pe ¯ L = . × ˆ n e ¯ n ¯ L . (21)To further simplify equation expressions, define α di = k B ¯ T i Ze ¯ B ¯ L ¯ V A , α de = k B ¯ T e e ¯ B ¯ L ¯ V A = α di Z τ i , (cid:15) p = k B ¯ T e m i ¯ V A , (cid:15) j = ¯ B µ e ¯ n ¯ L ¯ V A , (cid:15) v = m i Ze ¯ B ¯ t . (22)For simplicity, the over-hat symbol for normalized quantities is omitted for the rest of the paper, and the normalizedsix-field equations read ∂∂ t n i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ n i − n i B ˆ b × κ · ∇ ⊥ φ − α di B ˆ b × κ · ∇ ⊥ P i − n i B ∇ (cid:107) (cid:32) V (cid:107) i B (cid:33) + S n , (23) ∂∂ t (cid:36) = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ (cid:36) + B ∇ (cid:107) (cid:32) J (cid:107) B (cid:33) + (cid:15) p ˆ b × κ · ∇ P − (cid:15) p ˆ b × κ · ∇ π ci − α di B (cid:104) ˆ b × ∇ P i · ∇ (cid:16) ∇ ⊥ φ (cid:17)(cid:105) + (cid:16) ˆ b × ∇ n i · ∇ V E (cid:17) + α di B (cid:104) V E · ∇ (cid:16) ∇ ⊥ P i (cid:17) − ∇ ⊥ ( V E · ∇ P i ) (cid:105) , (24) ∂∂ t V (cid:107) i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ V (cid:107) i − (cid:15) p n i (cid:34) ∇ (cid:107) P + B ∇ (cid:107) (cid:32) π ci B (cid:33)(cid:35) − α di ˆ b n i B × ∇ P i · ∇ V (cid:107) i − V (cid:107) i S n n i , (25) ∂∂ t A ∗(cid:107) = −∇ (cid:107) φ + η (cid:107) ∇ ⊥ A (cid:107) + α de (cid:32) n e ∇ (cid:107) P e + . ∇ (cid:107) T e (cid:33) , (26) ∂∂ t T i = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) i ˆ b (cid:33) · ∇ T i − n i ∇ (cid:107) q (cid:107) i + (cid:15) q Z n i T / e (cid:32) T e τ i − T i (cid:33) − T i (cid:34)(cid:32) B ˆ b × κ (cid:33) · (cid:32) ∇ φ + α di ∇ P i n i + α di ∇ T i (cid:33) + B ∇ (cid:107) (cid:32) V (cid:107) i B (cid:33)(cid:35) − τ i π ci n i (cid:34) √ B ∇ (cid:107) (cid:16) √ BV (cid:107) i (cid:17) − α di n i B ˆ b · ∇ n i × ∇ T i (cid:35) − (cid:15) v B T i V (cid:107) i ˆ b × κ · ∇ V (cid:107) i + S Ei n i − T i S n n i , (27) ∂∂ t T e = − (cid:32) B ˆ b × ∇ ⊥ φ + V (cid:107) e ˆ b (cid:33) · ∇ T e − n e ∇ (cid:107) q (cid:107) e + .
71 23 (cid:15) j T e n e B ∇ (cid:107) (cid:32) J (cid:107) B (cid:33) − (cid:15) q Zn i T / e ( T e − τ i T i ) − T e (cid:34)(cid:32) B ˆ b × κ (cid:33) · (cid:32) ∇ φ − α de ∇ P e n e − α de ∇ T e (cid:33) + B ∇ (cid:107) (cid:32) V (cid:107) e B (cid:33)(cid:35) + (cid:15) p n e η (cid:107) J (cid:107) + S Ee n e − T e S n n e , (28)5ith (cid:36) = n i B (cid:32) ∇ ⊥ φ + n i ∇ ⊥ φ · ∇ ⊥ n i + α di ∇ ⊥ P i n i (cid:33) , (29) A ∗(cid:107) = (cid:16) − d e ∇ ⊥ (cid:17) A (cid:107) (30) π ci = − (cid:15) g T / i ˆ b × κ B · (cid:32) ∇ φ + α di ∇ P i n i (cid:33) + √ B ∂ (cid:107) ( √ BV (cid:107) i ) , (31) J (cid:107) = J (cid:107) + j (cid:107) where j (cid:107) = −∇ ⊥ A (cid:107) , (32)and V (cid:107) e = V (cid:107) i − (cid:15) j J (cid:107) n e . (33)In addition, to improve the numerical stability artificial forth order hyper-di ff usion terms are optional on the right-hand-side of the equations for each quantity. Depending on the nature of the problem and the resolution, these termscould be individually turned on to dissipate grid-scale turbulence eddies if necessary.
3. Coordinate, mesh, operators and boundary conditions
We briefly outline the coordinate system, mesh, operators and boundary conditions that are employed in theLandau fluid model in this section to lay the foundation for further discussion.
In the magnetic fusion community, there are more than a dozen choices of coordinate systems tailored for di ff erentresearch problems [21]; while in BOUT ++ , there are three coordinate systems are involved at di ff erent stages, namelythe 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 plasmaequilibrium; the flux coordinate system ( ψ, θ, ζ ) is widely used by theorists in analytical analysis; and the field-alignedcoordinate ( x , y , z ) is employed in the actual BOUT ++ simulation due to its superior numerical e ffi ciency. Becausethe cylindrical coordinate system ( R , φ, Z ) is mainly involved in the equilibrium calculation and mesh generation (e.g.,tokamak grid generator HYPNOTOAD [22]) 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 right-handed coordinate system as the one defined in EFIT [23] and kinetic EFIT. ( ψ, θ, ζ ) and ( x , y , z ), however, couldbe either right-handed or left-handed, depending on the poloidal field B θ direction. The positive toroidal ζ directionis defined as counter-clockwise when viewing from the top of torus, i.e., the opposite of the “normal” directionconsidered by experimentalists. The positive poloidal θ direction is defined as the direction of the polodial magneticfield B θ produced by a positive toroidal plasma current I p . In experiments, plasma current I p can be negative orpositive, 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 | e ψ | = R | B θ | , | e θ | = h θ , | e ζ | = R , (34)and hence | e ψ | = R | B θ | , | e θ | = h θ , | e ζ | = R . (35)The total magnetic field in a tokamak is composed by the poloidal and toroidal magnetic fields; conveniently, itcan be written as B = B θ + B ζ = ∇ ζ × ∇ ψ + F ( ψ ) ∇ ζ (36)with ψ the poloidal flux and F = RB ζ . The Jacobian of flux coordinate ( ψ, θ, ζ ) is thus simply J − = ∇ ψ × ∇ θ · ∇ ζ = ∇ ζ × ∇ ψ · ∇ θ = B θ · ∇ θ = B θ / h θ (37)6hich again can be either positive or negative depending on the sign of poloidal field σ B θ = B θ / | B θ | , or the directionof plasma current I p .Although the orthogonal flux coordinate ( ψ, θ, ζ ) works well in analytical analysis, BOUT ++ employs a field-aligned coordinate ( x , y , z ) to take the advantage of the strongly antiseptic feature of magnetized plasma by minimizingthe resolution requirement in the field-line direction: x = σ B θ ( ψ − ψ ) , y = θ, z = σ B θ (cid:32) ζ − (cid:90) θθ ν ( ψ, θ ) d θ (cid:33) . (38)Here ψ , θ are the reference poloidal flux and poloidal angle, normally set to be at the seperatrix and outboardmidplane, respectively. The local pitch of the magnetic field ν is defined as ν = B ζ h θ B θ R = ( F / R ) h θ B θ R = F J / R . (39)With Equation (38), one can write down the transformation of the contravariant and covaraint basis ( e i , e i ) betweentwo coordinate systems, e x e y e z = σ B θ − σ B θ I − σ B θ ν σ B θ e ψ e θ e ζ , e x e y e z = σ B θ σ B θ I ν σ B θ e ψ e θ e ζ , (40)where the integrated local shear I ( ψ, θ ) is defined as I ( ψ, θ ) = (cid:90) θθ ∂ν ( ψ, θ ) ∂ψ d θ. (41)As ∇ x = σ B θ ∇ ψ, ∇ y = ∇ θ, ∇ z = σ B θ ( ∇ ζ − I ∇ ψ − ν ∇ θ ) , (42)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 ) − , g i j = ( RB θ ) − I ( RB θ ) / h θ − σ B θ ν/ h θ − I ( RB θ ) − σ B θ ν/ h θ I ( RB θ ) + B / ( RB θ ) , g i j = RB θ ) + I R σ B θ B ζ h θ IRB θ IR σ B θ B ζ h θ IRB θ B h θ B θ σ B θ B ζ h θ RB θ IR σ B θ B ζ h θ RB θ R . (43)Similarly, one can write down the cross products of contravariant basis of field-aligned coordinate system, ∇ x × ∇ y = ( R | B θ | / h θ ) ˆ e ζ , ∇ x × ∇ z = − B , ∇ y × ∇ z = / ( Rh θ ) ˆ e ψ + IR | B θ | / h θ ˆ e ζ . (44)Note that in the field-aligned coordinate system, constant x and z means staying on same magnetic line, as themagnetic field could simply be represented in the Clebsch form B = ∇ z × ∇ x = B θ h θ (cid:16) e θ + ν e ζ (cid:17) = B θ h θ e y . (45)The above expression also indicates that the unit vector ˆ b = B / | B | = B θ / ( h θ | B | ) e y could be in either the same or theopposite direction of e y depending on the sign of B θ . The Jacobian of ( x , y , z ) is given by J − = ∇ x × ∇ y · ∇ z = B θ / h θ , (46)which is equivalent to Equation (37).As mentioned previously there are two types of magnetic field lines in the tokamak edge region. Figure 1 illustratesthe 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 magnetic7 igure 1: (a) Flux surfaces in 2D poloidal plane, and (b) magnetic field-lines in the 3D torus. White dashed lines denote separatrix, red, blue andgreen lines represents magnetic field-lines (b) and their poloidal projections (a) in the closed flux region, scrape-o ff -layer and private flux regionrespectively. field-lines on the poloidal plane form the so-called closed flux surfaces (red lines); while outside the separatrix, theprojections yield open flux surfaces as the magnetic field-lines are terminated at the divertor target plates. Dependingon whether the magnetic field-lines pass through the top poloidal plane, the open flux region is further categorizedinto the scrape-o ff -layer (blue lines) and the private flux region (green lines). With these three somewhat detachedregions, generating a mesh that could be parallelized e ffi ciently is no longer straightforward. However, by carefullydividing 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 capturethe geometric information, BOUT ++ ’s tokamak meshes are normally based on ideal MHD equilibrium either attainedfrom external solvers such as CORSICA or GATO, or reconstructed from experimental measurements via EFIT orkinetic EFIT. Auxiliary mesh generation codes (e.g., HYPERNOTOAD [22]) are developed to produce high accuracyBOUT ++ meshes with desired resolution and domain. Therefore, BOUT ++ is capable to simulate realistic tokamakdischarges 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, thegeneric representation of magnetic field configuration in BOUT ++ makes the Landau fluid model a flexible tool tostudy plasmas in a variety kinds of magnetic configurations, for example, the linear devices [24], the helical plasmasand even the stellarator plasmas (B Shanahan, B Dudson and P Hill, Plasma Phys. Control. Fusion 61 (2019) 025007). Here we list the most common operations involved in the fluid simulations generalized in the field-aligned coordi-nate system. In BOUT ++ , a vector v is by default in covariant form unless it is declared otherwise. However, it maybe transformed in contravariant form during certain operations to reduce the overall computation. The dot and cross products of two vectors u and v can be calculated either in covariant or contravariant forms, u · v = (cid:88) i , j g i j u i v j = (cid:88) i , j g i j u i v j , (47) u × v = J (cid:88) k ( u i v j − u j v i ) e k = J (cid:88) k ( u i v j − u j v i ) e k . (48)8 igure 2: Typical poloidal planes for (a) limiter, (b) single null and (c) double null divertor configurations and their corresponding xy meshes in theindex space. White dashed lines denote last-closed flux surface, or separatrix. Red, blue (and yellow), and green areas represent closed flux region,scrape-o ff -layer and private flux region. Gray dash lines indicate the locations of branch-cut. The gradient of a scalar f is most convenient to be evaluated in covariant components form ( u i ) as u = ∇ f = (cid:88) i ∂ f ∂ i e i . (49)In some cases, only the perpendicular component of ∇ f is needed which is defined as u ⊥ = ∇ ⊥ f = − ˆ b × ( ˆ b × ∇ f ) = − B B × ( B × ∇ f ) . (50)The divergence of a vector u is most convenient to be evaluated in contravariant components form ( u i ) as f = ∇ · u = J (cid:88) i ∂∂ i ( Ju i ) . (51)9he curl of a vector v is again most convenient to be evaluated in covariant components form ( v i ) and the result is incontravariant components form ( u i ), u = ∇ × v = J (cid:88) k (cid:32) ∂ v j ∂ i − ∂ v i ∂ j (cid:33) e k = J (cid:32) ∂ v z ∂ y − ∂ v y ∂ z (cid:33) e x + J (cid:32) ∂ v x ∂ z − ∂ v z ∂ x (cid:33) e y + J (cid:32) ∂ v y ∂ x − ∂ v x ∂ y (cid:33) e z . (52)The Laplacian operator, by definition, is ∇ f = J (cid:88) i ∂∂ i J (cid:88) j ∂ f ∂ j e j · e i , (53)or, in explicit form ∇ f = |∇ x | ∂ f ∂ x + |∇ y | ∂ f ∂ y + |∇ z | ∂ f ∂ z + ∇ x · ∇ y ∂ f ∂ x ∂ y + ∇ x · ∇ z ∂ f ∂ x ∂ z + ∇ y · ∇ z ∂ f ∂ y ∂ z + ∇ x ∂ f ∂ x + ∇ y ∂ f ∂ y + ∇ z ∂ f ∂ z . (54) In the fluid model, an important type of operator is called “Poisson bracket” defined as[ f , g ] = ˆ b B × ∇ f · ∇ g = ˆ b B · ∇ f × ∇ g , (55)which physically means convection processes (e.g., E × B convection if f = φ , or magnetic flutter e ff ects if f = A (cid:107) ).Rewriting ∇ f × ∇ g as ∇ f × ∇ g = (cid:32) ∂ f ∂ x ∂ g ∂ y − ∂ f ∂ y ∂ g ∂ x (cid:33) e x × e y + (cid:32) ∂ f ∂ x ∂ g ∂ z − ∂ f ∂ z ∂ g ∂ x (cid:33) e x × e z + (cid:32) ∂ f ∂ y ∂ g ∂ z − ∂ f ∂ z ∂ g ∂ y (cid:33) e y × e z , (56)and with the aid of Equations (44) and (45), the full expression for Poisson bracket becomes B × ∇ f · ∇ g = B ζ R | B θ | h θ (cid:32) ∂ f ∂ x ∂ g ∂ y − ∂ f ∂ y ∂ g ∂ x (cid:33) − B (cid:32) ∂ f ∂ x ∂ g ∂ z − ∂ f ∂ z ∂ g ∂ x (cid:33) + IB ζ R | B θ | h θ (cid:32) ∂ f ∂ y ∂ g ∂ z − ∂ f ∂ z ∂ g ∂ y (cid:33) . (57)Because E × B drift v E = E × B / B = B × ∇ φ/ B , E × B advection then becomes v E · ∇ g = B B × ∇ φ · ∇ g = [ φ, g ] . (58)In general, an advection process of quantity f with velocity u is easier to be evaluated in contravariant form as u · ∇ f = u x ∂ f ∂ x |∇ x | + u y ∂ f ∂ y |∇ y | + u z ∂ f ∂ z |∇ z | + (cid:32) u x ∂ f ∂ y + u y ∂ f ∂ x (cid:33) ∇ x · ∇ y (cid:32) u x ∂ f ∂ z + u z ∂ f ∂ x (cid:33) ∇ x · ∇ z + (cid:32) u y ∂ f ∂ z + u z ∂ f ∂ y (cid:33) ∇ y · ∇ z (in covariant form) = u x ∂ f ∂ x + u y ∂ f ∂ y + u z ∂ f ∂ z (in contravariant form) (59)The parallel derivative along an unperturbed magnetic field ˆ b = B / B is ∇ (0) (cid:107) f = ˆ b · ∇ f = B ( ∇ z × ∇ x ) · (cid:32) ∇ y ∂ f ∂ y (cid:33) = JB ∂ f ∂ y = B θ h θ B ∂ f ∂ y , (60)10nd the double derivative then becomes ∇ (cid:107) f = B θ h θ B ∂∂ y (cid:32) B θ h θ B ∂ f ∂ y (cid:33) = B θ h θ B ∂∂ y (cid:32) B θ h θ B (cid:33) ∂ f ∂ y + B θ h θ B ∂ f ∂ y . (61)The derivative along the normalized perturbed magnetic field ˜ b = ˜ B / B = (1 / B ) ∇ A (cid:107) × ˆ b is˜ b · ∇ f = (1 / B ) ∇ A (cid:107) × B · ∇ f = − (1 / B ) B × ∇ A (cid:107) · ∇ f = − [ A (cid:107) , f ] . (62)The last type of important operators in the toroidal system is the curvature operator ˆ b × ˆ κ · ∇ f , where ˆ κ is themagnetic line curvature, defined as ˆ κ = ˆ b · ∇ ˆ b = − ˆ b × ( ∇ × ˆ b ). The curvature operator could be rewritten in a similarform of the advection operator ˆ b × ˆ κ · ∇ f = v x ∂ f ∂ x + v y ∂ f ∂ y + v z ∂ f ∂ z (63)where the velocity components v i = ˆ b × ˆ κ · ∇ i with i = x , y , z are calculated during mesh generation and then passedalong as the prescribed inputs. Neglecting the perturbed magnetic field contribution is justified because normally it isthree 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, radial derivative inBOUT ++ in fact is often not done in field-aligned coordinate but “shifted” [13] and evaluated in a so-called “quasi-ballooning” coordinate. One may refer to References [1, 25] for details. 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 availablewithin BOUT ++ for di ff erent physics models, including the most commonly used Dirichlet, Neumann, Robin bound-ary conditions as well as the “relaxed” boundary which evolves the quantity towards the given boundary conditions atcertain rate in order to prevents transient discontinuities occurring at the boundaries. In general, for the global Landaufluid model, if applicable, Dirichlet boundary is applied to the perturbed quantities, and Neumann boundary is appliedto the equilibrium quantities to minimize the fluctuations at the boundary. The parallel y direction has two distinctboundary conditions depending on its corresponding magnetic topology. As illustrated in Figure 1, inside the sepa-ratrix magnetic field-lines wind around the major axis while simultaneously swirl around minor (magnetic) axis, andnest so-called closed flux surfaces. In most cases, a magnetic field-line will not end up at the same toroidal locationafter finishing one poloidal turn; therefore, a twist-shift operation [12, 13] is used to ensure the parallel continuity inthe closed flux region. In the open field-line region, i.e., the scrape-o ff -layer and the private flux region classical Bohmsheath criteria [26] are enforced at the boundary where magnetic field-lines are intercepted by limiter or divertor targetplates, V (cid:107) i = ± max[ C s , | V (cid:107) i | ] , (64) J (cid:107) = ± n e e (cid:34) C s − v th , e √ π exp (cid:32) − e Φ k B T e (cid:33)(cid:35) , (65) q (cid:107) e , i = ± γ e , i n e , i k B v the , i T e , i , (66)Here C s is the local sound speed, arising from the presheath electrostatic potential. In divertor configuration, theoryand model predict that ion could be supersonic near the divertor plates due to the Laval nozzle e ff ect [27, 28], su-personic ion flow is also occasionally observed in UEDGE and BOUT ++ transport simulations [29]. Therefore, ionparallel velocity V (cid:107) i is compared to the the local sound speed C s at each time step – if it is less than C s , then | V (cid:107) 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 coe ffi cient, local electron and ion temperature ratio, target plate floating potential [30].The typical values used in BOUT ++ simulations are γ e (cid:39) , γ i (cid:39) .
5. As will be discussed in next section, theimplementation of this boundary condition depends on the specific parallel heat flux model to be used. In addition the11ollowing boundary conditions on density, vorticity and electrostatic potential are applied to eliminate the potentialbuild-up of large discontinuities at the boundary, ∇ (cid:107) n i = ∇ (cid:107) (cid:36) = , (67) ∇ (cid:107) φ = − η (cid:107) j (cid:107) + ∇ (cid:107) P e en e + . k B e ∇ (cid:107) T e . (68)
4. Parallel heat flux models
The Braginskii transport model [15] derived under the assumption that plasma is magnetized and collsional re-quires that L ⊥ (cid:29) ρ and L (cid:107) (cid:29) λ mfp . Here ρ is the particle (electron and ion) gyroradius, λ m f p is the particle mean freepath, and L ⊥ , L (cid:107) are the characteristic lengths in perpendicular and parallel directions. Under the strong collisionallimit, the well-known Braginskii (or, Spitzer-H¨arm) heat flux model is deduced by assuming the perturbed particledistribution slightly deviates from the equilibrium Maxwellian distribution, q B (cid:107) j = − κ j (cid:107) ∇ (cid:107) (cid:16) k B T j (cid:17) with κ e (cid:107) = . n e k B T e τ ei m e , κ i (cid:107) = . n i k B T i τ ii m i (69)which implies that heat flux depends on local temperature gradient, and heat (or, energy) is transported solely throughcollisions (or, random walk process) among particles.Recent particle simulations indicate that Braginskii parallel heat flux model is reasonably accurate (i.e., less than10% error) when the parallel inhomogeneity length is at least on the two order of magnitude longer than the particlemean free path (i.e., L (cid:107) /λ m f p ≥ (10 )) and becomes invalid when L (cid:107) /λ m f p ≤ (10) [31].In the tokamak edge region, L ⊥ is the equilibrium length such as L n , L T which are typically a few centimeters,while L (cid:107) ∼ π qR in large aspect ratio limit is on the order of tens meters where q ≈ ∼ R is the major radius of the machine. Assuming that plasma temperature is isotropic (i.e., T ⊥ = T (cid:107) = T ),then ρ = v th /ω c ∝ m / T / and λ m f p ∼ v th τ ∝ T / n is independent of mass. As modern and future tokamaksproduce higher and higher temperature plasmas, the latter requirement ( L (cid:107) (cid:29) λ mfp ) is no long be satisfied due tothe strong temperature dependence of λ m f p . For example, Table 1 examines the validity of Braginskii model for theH-mode edge of three tokamaks: Alcator C-Mod (compact high field, relatively high collisionality machine), DIII-D(medium size, relatively low collisionality machine) and ITER (future experimental reactor). These estimate indicatethat Braginksii parallel heat flux model merely holds near the separatrix and breaks down at the pedestal region. Thedirect consequence of violation of strong collisionality criterion L (cid:107) (cid:29) λ mfp is that the classical Braginskii (or, Spitzer-H¨arm) parallel heat flux model (Equation 69) is no longer accurate and overestimates the parallel heat flux by ordersof magnitude in the weakly collisional regime as demonstrated in Figure 3. Although the classical Braginskii parallelheat flux model is still widely used for more collisional L-mode edge plasma simulations nowadays, [32, 33] a moreadvanced means of parallel heat flux evaluation for weakly collsional plasma is needed in fluid simulations. Intuitively, in a collisionless plasma, parallel heat transport is mainly done by particles freely moving along themagnetic field-lines which also sets a upper limit of heat flux. In fluid description, this free streaming heat flux isdefined as q FS (cid:107) = nk B T v th where v th is the thermal speed. Hence, a flux-limited expression for parallel heat flux isconstructed to constrain the unrealistically large value predicted by Braginskii model [34], q FL (cid:107) = q B (cid:107) + α q FS (cid:107) − (70)where α ∈ [0 . ,
1] is a parameter to be adjusted in di ff erent scenarios. However, Equation (70) becomes singular atwhere q B (cid:107) + α q FS (cid:107) ≈ q FS (cid:107) is always positive while q B (cid:107) can be negative depending on local temperaturegradient and | q B (cid:107) | (cid:29) q FS (cid:107) in the weakly collisional regime. In practice, q FS (cid:107) is replaced with the “di ff usive” free12 able 1: Braginskii validation for Deuterium plasma Mid-pedestal SeparatrixCMod DIII-D ITER CMod DIII-D ITER B ( T ) 5 2 5 5 2 5 q R ( m ) 0.68 1.67 6.2 0.68 1.67 6.2 n (10 m − ) 1.6 0.6 0.6 0.8 0.4 0.3 T ( eV ) 150 250 2000 50 80 200 L ⊥ ( cm ) 1 1.5 3 2 3 6 L (cid:107) ( m ) 17 42 117 21 52 156 ρ e ( cm ) 5 . × − . × − . × − . × − . × − . × − ρ i ( cm ) 3 . × − . × − . × − . × − . × − . × − λ m f p ( m ) 1.4 9.5 505 0.34 1.6 12.1streaming heat flux q FS(d) (cid:107) = − nv th L (cid:107) ∇ (cid:107) ( k B T ) which ensures it has the same sign as q B (cid:107) at all locations. Consequently,the resulting flux-limited parallel heat flux is again dependent on local temperature gradient, similar to the Braginskiimodel. It can be rewritten as q FL (cid:107) j = − κ j e ff ∇ (cid:107) ( k B T j ) with the e ff ective thermal conductivity κ j e ff = α n j v th , j L (cid:107) κ j (cid:107) α n j v th , j L (cid:107) + κ j (cid:107) . (71) Although the flux-limited parallel heat flux model overcomes the overestimation of Braginskii model, it is stillan ad hoc model. A more rigorous approach is to develop a theory-based parallel heat flux model, for example, theLandau fluid model for weakly collisional plasmas [35] q LF (cid:107) j , k = − n j k B (cid:114) π v th , j ik (cid:107) | k (cid:107) | + α j /λ j T j , k , (72)where α e (cid:39) .
499 and α i (cid:39) . ff erent coe ffi cients and an extension of Hammett-Perkins closure [38] – thefirst and the simplest Landau fluid closure that describes plasma kinetic responses in the collsionless limit, q HP (cid:107) , k = − n j k B (cid:114) π v th , j ik (cid:107) | k (cid:107) | T j , k (73)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 e ff ects [39], this closure predicts aslightly lower heat flux (with a maximum di ff erence by 25%) in the intermediate collisionality regime.Implementing such a closure in BOUT ++ turbulence model is numerically challenging. This is because Landaufluid closures are usually entail wave-vector space (i.e., k − space) and implicitly assume periodicity along magneticfield-lines; while BOUT ++ turbulence model calculates quantities in configuration space. A novel solution to resolvethese issues is the fast non-Fourier method [40] which approximates parallel heat flux with summation of Lorentiziansas q (cid:107) j , k ≈ N − (cid:88) n = q ( n ) (cid:107) j , (74)where q ( n ) (cid:107) j = − n j k B (cid:114) π v th , j α n ξ j k (cid:107) ( ξ j k (cid:107) ) + β n iT j , k . (75)13 igure 3: Parallel heat flux v.s. collisionality. Here ξ j = λ m f p , j /α j and ( α n , β n ) are optimized fitting coe ffi cients of xx + ≈ N − (cid:88) α n xx + β n (76)as the ones given in Table 2 of Reference [41]. It is therefore straightforward to transform the above Lorentizians backto configuration space by replacing ik (cid:107) and k (cid:107) with ∇ (cid:107) and −∇ (cid:107) , (cid:16) β n − ξ j ∇ (cid:107) (cid:17) q ( n ) (cid:107) j = − n j k B (cid:114) π v th , j ξ j α n ∇ (cid:107) T j , (77)and the evaluation of parallel heat flux in k − space becomes solving a series of boundary value problems in configura-tion space. The fast non-Fourier method has been developed and validated with a one-dimensional uniform grid [41].As illustrated in Figure 4, it attains more accurate approximation and reaches further into the collisionless regime byadding more Lorentizians. In most tokamak edge simulations, N = ffi cient as it is valid for ξ k (cid:107) ∼ k (cid:107) λ m f p < with a relative error about 2%. One can further push the applicable limit to k (cid:107) λ m f p < and / or reduce the relativeerror to about 1% by using N =
12 Lorentizians. This treatment of Landau fluid closure has now been extended andimplemented 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 Landaufluid 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 bothscenarios, electron parallel heat flux is evaluated assuming there is a relatively large (30%) temperature perturbationis 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 regionwhere temperature gradient is finite (thus, “local” model) and its amplitude is about two orders of magnitude largerthan either flux-limited or the Landau fluid estimate. The flux-limited model seems able to e ff ectively scale down theheat 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”e ff ect as the predicted parallel heat flux extends to the region where local temperature gradient vanishes (e.g., thetokamak bottom region). More interestingly, our tests also demonstrate the transition from local-dominant transportto nonlocal-dominant transport. In C-Mod plasma, parallel heat flux still marginally depends on the local temperaturegradient (Figures 5(f)). As the plasma becomes less collisional, this local e ff ect weakens and the nonlocal e ff ectdominates 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 am-bipolar sheath heat flux boundary condition q sheath (cid:107) at the divertor targets (i.e., Equation 66) as shown in Figures 5(i)14 igure 4: Accuracy and the validity of the Landau fluid closure v.s. collisionality for di ff erent numbers of Lorentizians. and 6(i). This is again due to the distinct inherent transport characteristics. Because Braginskii and flux-limited mod-els are local, q sheath (cid:107) is actually not enforced on parallel heat flux but on the plasma temperature as ∇ (cid:107) T = ± q sheath (cid:107) /κ which tends to produce a steep parallel gradient at the divertor targets. While for Landau fluid closure, q sheath (cid:107) isprovided 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 paralleldynamics may also be altered significantly. Perhaps a more important quantity is the parallel gradient of parallel heatflux ( ∇ (cid:107) q (cid:107) ) the since it directly influences temperature evolution. Our initial analysis suggests that Landau fluid closureresults in a nearly one order of magnitude weaker heat conduction comparing to the value given by the flux-limitedmodel both in the closed flux region and in the SOL, as attested by the substantial diminishing amplitudes of ∇ (cid:107) q (cid:107) (reddotted 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 [42] doesn’t fit in the scope of this paper and will be left for afuture study.Although the Landau fluid closure is arguably the most sophisticated treatment of parallel heat flux in nowadaysfluid edge turbulence models, it is worth to clarify that there are still a few subtleties in this approach. For instance,the temperature anisotropy e ff ect is neglected in order to be consistent with Braginskii transport equations, and theclassical Bohm’s sheath criteria for collisional plasma is used at the boundaries in the open field-line region. Theformer can be addressed in a more complete BOUT ++ gyro-Landau-fluid turbulence model. Although one may arguethat 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 = √ T / m i π vs c s ) because ions that reach the divertorplate will escape. While the heat flux predicted by Landau fluid closure (72) is found to have a good agreement withthe result from a Fokker-Planck code with a small amplitude perturbation in weakly collisional regime [35], a furthervalidation with extended parameter regime will be performed and reported in future publications. For the sake of completeness, the normalization of di ff erent parallel heat flux models are listed below. In BOUT ++ transport and turbulence codes, q (cid:107) j is naturally normalized with ¯ q j = ¯ nk B ¯ V A ¯ T j and the normalized parallel heat fluxfor Braginskii expression and flux-limited model along with their thermal conductivities are q B (cid:107) j = − κ j (cid:107) ∇ (cid:107) T j and q FL (cid:107) j = − κ j e ff ∇ (cid:107) T j (78)with κ e (cid:107) = . τ ei T / e v th , e , κ e (cid:107) = . τ ii T / i v th , i , κ j e ff = α n j v th , j L (cid:107) κ j (cid:107) α n j v th , j L (cid:107) + κ j (cid:107) . (79)15 igure 5: Comparison of parallel heat flux for typical C-Mod H-mode parameters: singly charged deuterium plasma equilibrium (a) density n , (b)electron temperature T e , (c) perturbed electron temperature T e and the parallel heat flux predictions from (d) Braginskii model, (e) flux-limitedmodel ( α = . Here the plasma density and temperature n j , T j , thermal speed v th , j and parallel characteristic length L (cid:107) are normalizedto ¯ n , ¯ T j , ¯ V A and ¯ L respectively.For the Landau fluid closure, Equation (77) that solves Lorentizians of parallel heat flux after normalizationbecomes (cid:16) β n − ξ j ∇ (cid:107) (cid:17) q ( n ) (cid:107) j = − n j (cid:114) π v th , j ξ j α n ∇ (cid:107) T j , (80)where n j , T j , ξ j and v th , j are again normalized to ¯ n , ¯ T j , ¯ L and ¯ V A .
5. Field solver
In the drift-reduced fluid formalism, Laplacian inversion is required to evaluate electrostatic potential φ fromgeneral vorticity (cid:36) , and parallel magnetic potential A (cid:107) (or perturbed magnetic flux ψ = − A (cid:107) ) from modified magneticflux A ∗(cid:107) if electron inertia is included in the model, ∇ ⊥ φ + n i ∇ ⊥ φ · ∇ ⊥ n i = B (cid:36) n i − α di n i ∇ ⊥ P i , (81) d e ∇ ⊥ A (cid:107) − A (cid:107) = − A ∗(cid:107) . (82)Therefore, Laplacian inversion solvers are designed in BOUT ++ to solve both equations in a generalized form d ∇ ⊥ f + c ∇ ⊥ c · ∇ ⊥ f + a f = b . (83)16 igure 6: Comparison of parallel heat flux for typical DIII-D H-mode parameters: singly charged deuterium plasma equilibrium (a) density n , (b)electron temperature T e , (c) perturbed electron temperature T e and the parallel heat flux predictions from (d) Braginskii model, (e) flux-limitedmodel ( α = . Since BOUT ++ employs the field aligned coordinate system, the perpendicular component of Laplacian operatorused in the vorticity equation and Ohm’s law, by definition, is ∇ ⊥ f = ( ∇ − ∇ (cid:107) ) f (84) = G x ∂ f ∂ x + G y ∂ f ∂ y + G z ∂ f ∂ z + g xx ∂ f ∂ x + g yy ∂ f ∂ y + g zz ∂ f ∂ z + g xy ∂ f ∂ x ∂ y + g xz ∂ f ∂ x ∂ z + g yz ∂ f ∂ y ∂ z − g yy ∂ f ∂ y − J ∂∂ y (cid:32) Jg yy (cid:33) ∂ f ∂ y where G i = ∇ i for i = x , y , z and the ∂ xy term is omitted as g xy = ++ separates theunknown quantify f into non-axisymmetric ( n (cid:44)
0) and axisymmetric ( n =
0) components to improve the numericalperformance by reducing the dimensionality of the problem. Here n represents the toroidal mode number. (cid:44) ) field For the non-axisymmetric field, considering the strongly anisotropic nature of magnetized plasmas, the parallelderivatives are normally several orders of magnitude smaller than the perpendicular derivatives, i.e., k (cid:107) (cid:28) k ⊥ . Thereforewe may neglect any ∂/∂ y terms and also the G x , G z terms (as the grid is normally uniform in x and z direction), and17 igure 7: q (cid:107) and ∇ (cid:107) q (cid:107) predicted by three models along field-lines (a) at pedestal top ( ψ N = .
95) and (b) in the SOL ( ψ N = .
01) in the DIII-D test(i.e., Figure 6).
Equation (84) becomes ∇ ⊥ f = g xx ∂ f ∂ x + g zz ∂ f ∂ z + g xz ∂ f ∂ x ∂ z = ( RB θ ) (cid:40) ∂ f ∂ x − I ∂ f ∂ x ∂ z + (cid:34) I + B ( RB θ ) (cid:35) ∂ f ∂ z (cid:41) . (85)Rewriting above equation in flux coordinate system to drop all the I terms and then applying Fourier transformalong the toroidal direction, yields a 1D boundary value problem and can be solved e ffi ciently with the tridiagonalsolvers [25]. = ) field However, the above approach doesn’t apply to the axisymmetric k z = n =
0) zonal field component. Intu-itively, this means that the compressed 1D profile ( x distribution) does not contain enough information to reconstructa 2D structure ( xy − plane); therefore, one needs to include y variation in order to correctly solving the zonal field com-ponent. In other similar edge turbulence codes, this process is normally done with a hard-coded algebraic multigridsolver [43, 44, 45]. BOUT ++ utilizes external libraries (e.g., PETSc [46], HYPRE [47]) to solve this problem thus hasa more flexible option on the numerical algorithm to be used [48, 25]. In BOUT ++
2D transport code, the relaxationmethod is used to solve the zonal potential equation for a steady-state solution. [49]For n = ∂/∂ z terms vanishes, and Equation (84) becomes ∇ ⊥ f = G x ∂ f ∂ x + G y ∂ f ∂ y + g xx ∂ f ∂ x + g yy ∂ f ∂ y − g yy ∂ f ∂ y − J ∂∂ y (cid:32) Jg yy (cid:33) ∂ f ∂ y . (86)Above equation is then discretized either based on the second order finite di ff erence scheme (thus, 5-point stencil) orthe forth order finite di ff erence 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 per-formed. That is, giving an self-consistent set of profiles ( n i , T i , B ) and assuming (cid:36) =
0, evaluate the right-hand-sideof Equation (81), i.e., ∇ ⊥ φ + n i ∇ ⊥ φ · ∇ ⊥ n i = − α di n i ∇ ⊥ P i = RHS ; (87)and calculate φ by solving the 2D boundary value problem; then re-evaluate the left-hand-side of Equation (81) (de-noted as RHS (cid:48) ) based on φ ; and finally the absolute error δ RHS = RHS (cid:48) − RHS is used to quantify the performanceof the solver in this quasi-realistic setting. In this test, Equation (81) is discretized with a 5-point stencil, and thegeneralized minimal residual (GMRES) method and the successive over-relaxation (SOR) method in the PETSc li-brary are picked as the matrix solver and the preconditioner. Dirichlet boundary condition ( φ = constant) is appliedon the outside boundary (wall-side) and Neumann boundary condition ( ∂φ/∂ψ = ∝ V E , p ) is enforced on the insideboundary (core-side) to eliminate any external poloidal angular momentum sourcing. Additional options, such as the18bsolute and relative error tolerances (cid:15) A , (cid:15) 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 , δ RHS ∼ − (Figure 8(f)) is about eight orders of magnitude smaller than the initial input (Figure 8(c)), suggestingthe the axisymmetric field solver is correctly implemented. Moreover, the resulting error δ RHS tends to collocate atwhere
RHS is finite.A convergence test is carried out by varying the relative error tolerance (cid:15) R with and without the ∇ ⊥ n i · ∇ ⊥ φ crossterm on the left-hand-side of Equation (81), and the result is summarized in Figure 10(a). Not surprisingly, theaxisymmetric field solver converges. As (cid:15) R decreases, the maximum absolute error of the test | δ RHS | max decreaseslinearly as well till it saturates at 10 − level. As expected, the solver slows down as (cid:15) R becomes smaller, i.e., thenumber of iterations required to achieve certain criterion increases. However, the apparent discrepancy between thetwo red lines indicates that the number of iterations is reduced substantially when the cross term is neglected (i.e.,applying Boussinesq approximation – a widely used but not fully justified assumption in reduced-fluid research)especially with a fixed relatively small (cid:15) R . It is worth to point out that the exact number of iterations depends on quitea few factors, for example, grid resolution, smoothness and consistency of the profiles ( RHS and n i ), choices of matrixsolver 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 , φ and reproduce the right-hand-side of Equation (81) in a good precisionalthough in this case, error seems to concentrate at the far private flux region (likely due to the choices of initialprofiles and boundary conditions). Interestingly, in the convergence test, no significant discrepancy of the number ofiterations with and without the cross term is observed. Figure 8: Inverse-and-forward electrostatic potential test of the axisymmetric solver in a shift-circular configuration: (a) density, (b) ion temperatureand (c) evaluated
RHS profiles, (d) calculated φ , (e) re-evaluated RHS (cid:48) and (f) the absolute error δ RHS . This test is done with cross term and (cid:15) R = − . The new axisymmetric field solver is further benchmarked against relaxation method [49] currently implementedin BOUT ++ ’s transport model. As plotted in Figure 11, for a given initial potential and vorticity profiles φ , (cid:36) ,both methods are able to recover the initial potential profile. The overall small discrepancy between two solutions islikely due to the slightly shifted location of boundary condition implementation – while relaxation method enforcesboundary condition at the boundary grid points, the axisymmetric field solver enforces boundary condition at theboundary cell centers.Authors would like to point out that the new field solver is able to invert the a more accurate form of vorticity,19 igure 9: Inverse-and-forward electrostatic potential test of the axisymmetric solver in a C-Mod-like lower single null divertor configuration: (a)density, (b) ion temperature and (c) evaluated RHS profiles, (d) calculated φ , (e) re-evaluated RHS (cid:48) and (f) the absolute error δ RHS . This test isdone with cross term and (cid:15) R = − .Figure 10: The maximum absolute error (blue lines) and the number of iterations (red lines) in the inverse-and-forward tests v.s. relative tolerance (cid:15) R with (dashed lines) and without (solid lines) cross term with (a) shift circular and (b) lower single null configuration. e.g., including the n = ∇ ⊥ φ · ∇ ⊥ n i and / or the iondensity dependence in front of ∇ ⊥ φ term (i.e., so called Boussinesq approximation) in order to simplify the expression20 igure 11: Initial (a) electrostatic potential and corresponding (b) vorticity v.s. the solutions from (c) direct axisymmetric field solver and (d)relaxation method for a lower single null configuration. of vorticity, hence the numerical calculation, (cid:36) ≈ n i B (cid:32) ∇ ⊥ φ + α di n i ∇ ⊥ P i (cid:33) (88) (cid:36) ≈ B (cid:16) ∇ ⊥ φ + α di ∇ ⊥ P i (cid:17) (89)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. Our preliminarytests indicate that the cross term and the 1 / n i dependence do influence the inverted axisymmetric electrostatic poten-tial and the corresponding radial electric field (i.e., zonal field) as shown in Figure 12. The impact of Boussinesqapproximation on edge turbulence and flow generation simulation is subject of future research.
6. Flux driven source
In the edge turbulence simulations, plasma is constantly transported outward once instability is excited, it istherefore reasonable to refuel particles and energy from the core side in order to prevent the collapse of densityand temperature profiles in the transport time scale simulations. Meanwhile, plasma sources are also necessary tosimulate certain transient events, e.g., the thermal quench phase of tokamak disruption when the hot core plasma ispushed towards edge in a very short time-scale ( 100 µ s ).21 igure 12: Solutions of axisymmetric electrostatic potentials and the corresponding radial electric fields assuming vorticity (cid:36) = ff erentvorticity expression: (a,d) complete form (29), (b,e) without cross term (88), (c,f) without density dependence (88). The newly added sourcing options therefore largely expands the range of physics problems that BOUT ++ edgemodel is able to investigate. A good demonstration of this new capability is the preliminary study of divertor heat loadresponse to excessive heating which is an important issue in tokamak disruption mitigation.A nonlinear electromagnetic simulation is first carried out in a DIII-D like H-mode lower single null plasmawithout sourcing for the first 0.14 ms . While for the next 80 µ s , 1 GW power, or 80 kJ energy in total, equallypartitioned between electrons and ions, is then injected to the pedestal top region (0 . < ψ N < .
96) to mimicintensive energy flux coming from the core in the thermal quench phase. This is equivalent to roughly 15% of thestored plasma thermal energy in a typical DIII-D H-mode discharge is deposited at the edge region within 80 µ s . At t = . ms , energy source is turned o ff . Figure 13 plots the time history of outboard midplane electron temperature T e , OMP at the pedestal top ( ψ N = .
95) and the corresponding peak magnitude of the heat flux at outer divertor targetplate. A monotonic, almost linear increasing of T e , OMP is observed once the power injection starts, followed by amuch slower decreasing phase once the source is turned o ff . The divertor heat flux, on the other hand, takes a while toresponse to the change of pedestal temperature. The 30 µ s lag (e.g., between the peak temperature and peak divertorheat load) is mainly due to the parallel conduction process where the plasma inside the separatrix first enters the SOLprimarily near the outboard midplane and then is transported to the divertor target. In this run, doubled pedestal toptemperature yields more than 40 times larger heat flux. The dramatic increasing of heat load is caused by the enhancedturbulence activity at the tokamak edge as illustrated in Figure 14. Before the intensive energy injection, the systemis marginally stable to the intermediate- n peeling-ballooning modes localized at the steep gradient region; once theheating is on, elevated temperature profile starts to excite instabilities and eventually drives whole system turbulent.Although the underlying physics in process is yet to be explored, the versatile flux-driven source option provides aunique opportunity to study edge plasma dynamics in this pulsed over-driven setting.22 igure 13: Outboard midplane electron temperature at ψ N = .
95 and maximum outer divertor heat load evolution.Figure 14: Poloidal snapshots of normalized plasma pressure perturbation at three di ff erent time, (a) t = . ms , i.e., the beginning of energyinjection, (b) t = . ms , i.e., the end of energy injection and (c) t = . ms which is roughly at when the peak heat load is reached.
7. Summary
The drift-reduced two fluid six-field turbulence model in BOUT ++ framework is perhaps the most widely usedtokamak edge turbulence code in the magnetic fusion community. It evolves plasma density, electron and ion temper-ature, ion parallel velocity, electrostatic potential and parallel magnetic flux in a realistic tokamak magnetic configura-tion. In this paper, we present in detail the current status of this model in a self-contained manner, including the modelequations, normalization, coordinate system, operators, boundary conditions and the recently added features such asflux-driven source, Landau fluid closure on parallel heat flux and the axisymmetric field solver. These new featuresgreatly 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.23t 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 underdevelopment; while on the physics side, a comprehensive fluid neutral model and a ML-enabled surrogate model [50]of an even more sophisticated kinetic Landau fluid closure are also on the way. Acknowledgments
We thank Drs. J. Chen, A. Dimits, B. Dudson, N. Li, C. Ma, J. Omotani and M. Umansky for helpful discussions.This work is performed under the auspices of the U.S. Deportment of Energy (US DOE) by Lawrence LivermoreNational Laboratory under Contract DE-AC52-07NA27344 and supported by the US DOE Tokamak Disruption Sim-ulation (TDS) Scientific Discovery Through Advanced Computing (SciDAC) program. This research used resourcesof the National Energy Research Scientific Computing Center, a DOE O ffi ce of Science User Facility supported bythe O ffi ce of Science of the US DOE under Contract No. DE-AC02-05CH11231. LLNL-JRNL-812150. References [1] B. Dudson, M. Umansky, X. Xu, P. Snyder, H. Wilson, Bout ++ : A framework for parallel plasma fluid simulations, Computer PhysicsCommunications 180 (9) (2009) 1467–1480.[2] P. Ricci, F. Halpern, S. Jolliet, J. Loizu, A. Mosetto, A. Fasoli, I. Furno, C. Theiler, Simulation of plasma turbulence in scrape-o ff layerconditions: the gbs code, simulation results and code validation, Plasma Physics and Controlled Fusion 54 (12) (2012) 124047.[3] P. Tamain, H. Bu ff erand, G. Ciraolo, C. Colin, D. Galassi, P. Ghendrih, F. Schwander, E. Serre, The tokam3x code for edge turbulence fluidsimulations of tokamak plasmas in versatile magnetic geometries, Journal of Computational Physics 321 (2016) 606–623.[4] B. Zhu, M. Francisquez, B. N. Rogers, Gdb: A global 3d two-fluid model of plasma turbulence and transport in the tokamak edge, ComputerPhysics Communications 232 (2018) 46–58.[5] A. Stegmeir, D. Coster, A. Ross, O. Maj, K. Lackner, E. Poli, Grillix: a 3d turbulence code based on the flux-coordinate independent approach,Plasma Physics and Controlled Fusion 60 (3) (2018) 035005.[6] B. Scott, Edge turbulence and its interaction with the equilibrium, Contributions to Plasma Physics 46 (7-9) (2006) 714–725.[7] C. Ma, X. Xu, Global kinetic ballooning mode simulations in bout ++ , Nuclear Fusion 57 (1) (2016) 016002.[8] C.-S. Chang, S. Ku, H. Weitzner, Numerical study of neoclassical plasma pedestal in a tokamak geometry, Physics of Plasmas 11 (5) (2004)2649–2667.[9] A. H. Hakim, N. R. Mandell, T. Bernard, M. Francisquez, G. Hammett, E. Shi, Continuum electromagnetic gyrokinetic simulations ofturbulence in the tokamak scrape-o ff layer and laboratory devices, Physics of Plasmas 27 (4) (2020) 042304.[10] M. Dorf, M. Dorr, Progress with the 5d full-f continuum gyrokinetic code cogent, Contributions to Plasma Physics (2020) e201900113.[11] X. Xu, R. H. Cohen, Scrape-o ff layer turbulence theory and simulations, Contributions to Plasma Physics 38 (1-2) (1998) 158–170.[12] X. Xu, M. Umansky, B. Dudson, P. Snyder, Boundary plasma turbulence simulations for tokamaks, Communications in ComputationalPhysics 4 (5) (2008) 949–979.[13] M. Umansky, X. Xu, B. Dudson, L. LoDestro, J. Myra, Status and verification of edge plasma turbulence code bout, Computer PhysicsCommunications 180 (6) (2009) 887–903.[14] T. Xia, X. Xu, P. Xi, Six-field two-fluid simulations of peeling–ballooning modes using bout ++ , Nuclear Fusion 53 (7) (2013) 073009.[15] S. Braginskii, Transport processes in a plasma, Reviews of plasma physics 1 (1965).[16] A. N. Simakov, P. J. Catto, Drift-ordered fluid equations for field-aligned modes in low- β collisional plasma with equilibrium pressurepedestals, Physics of Plasmas 10 (12) (2003) 4744–4757.[17] J. Wesson, D. J. Campbell, Tokamaks, Vol. 149, Oxford university press, 2011.[18] T. Xia, X. Xu, Five-field simulations of peeling-ballooning modes using bout ++ code, Physics of Plasmas 20 (5) (2013) 052102.[19] P. Xi, X. Xu, T. Xia, W. Nevins, S. Kim, Impact of a large density gradient on linear and nonlinear edge-localized mode simulations, NuclearFusion 53 (11) (2013) 113020.[20] Z. Wang, X. Xu, T. Xia, T. Rognlien, 2d simulations of transport dynamics during tokamak fuelling by supersonic molecular beam injection,Nuclear Fusion 54 (4) (2014) 043019.[21] O. Sauter, S. Y. Medvedev, Tokamak coordinate conventions: Cocos, Computer Physics Communications 184 (2) (2013) 293–302.[22] B. D. Dudson, A. Allen, G. Breyiannis, E. Brugger, J. Buchanan, L. Easy, S. Farley, I. Joseph, M. Kim, A. McGann, et al., Bout ++ : Recentand current developments, Journal of Plasma Physics 81 (1) (2015).[23] L. Lao, H. S. John, R. Stambaugh, A. Kellman, W. Pfei ff er, Reconstruction of current profile parameters and plasma shapes in tokamaks,Nuclear fusion 25 (11) (1985) 1611.[24] P. Popovich, M. Umansky, T. Carter, B. Friedman, Analysis of plasma instabilities and verification of the bout code for the large plasmadevice, Physics of Plasmas 17 (10) (2010) 102107.[25] H. Seto, X. Xu, B. D. Dudson, M. Yagi, Interplay between fluctuation driven toroidal axisymmetric flows and resistive ballooning modeturbulence, Physics of Plasmas 26 (5) (2019) 052507.[26] P. C. Stangeby, The plasma boundary of magnetic fusion devices, CRC Press, 2000.[27] R. H. Cohen, D. Ryutov, Drifts, boundary conditions and convection on open field lines, Physics of Plasmas 6 (5) (1999) 1995–2001.[28] S. Togo, T. Takizuka, D. Reiser, M. Sakamoto, N. Ezumi, Y. Ogawa, K. Nojiri, K. Ibano, Y. Li, Y. Nakashima, Self-consistent simulation ofsupersonic plasma flows in advanced divertors, Nuclear Fusion (2019).
29] T. Rognlien, N. Li, personal communication (2019).[30] P. Stangeby, Plasma sheath transmission factors for tokamak edge plasmas, The Physics of fluids 27 (3) (1984) 682–690.[31] Z. Guo, X.-Z. Tang, C. McDevitt, Parallel heat flux and flow acceleration in open field line plasmas with magnetic trapping, Physics ofPlasmas 21 (10) (2014) 102512.[32] B. Zhu, M. Francisquez, B. N. Rogers, Up–down symmetry breaking in global tokamak edge simulations, Nuclear Fusion 58 (10) (2018)106039.[33] W. Zholobenko, A. Stegmeir, T. Body, A. Ross, P. Manz, O. Maj, D. Coster, F. Jenko, M. Francisquez, B. Zhu, et al., Thermal dynamics inthe flux-coordinate independent turbulence code grillix, Contributions to Plasma Physics 60 (5-6) (2020) e201900131.[34] J. Omotani, B. Dudson, Non-local approach to kinetic e ff ects on parallel transport in fluid models of the scrape-o ff layer, Plasma Physics andControlled Fusion 55 (5) (2013) 055009.[35] M. Umansky, A. Dimits, I. Joseph, J. Omotani, T. Rognlien, Modeling of tokamak divertor plasma for weakly collisional parallel electrontransport, Journal of Nuclear Materials 463 (2015) 506–509.[36] M. A. Beer, G. W. Hammett, Toroidal gyrofluid equations for simulations of tokamak turbulence, Physics of Plasmas 3 (11) (1996) 4046–4064.[37] P. Snyder, G. Hammett, W. Dorland, Landau fluid models of collisionless magnetohydrodynamics, Physics of Plasmas 4 (11) (1997) 3974–3985.[38] G. W. Hammett, F. W. Perkins, Fluid moment models for landau damping with application to the ion-temperature-gradient instability, Physicalreview letters 64 (25) (1990) 3019.[39] L. Wang, X. Xu, B. Zhu, C. Ma, Y.-a. Lei, Deep learning surrogate model for kinetic landau-fluid closure with collision, AIP Advances 10 (7)(2020) 075108.[40] A. Dimits, I. Joseph, M. Umansky, A fast non-fourier method for landau-fluid operators, Physics of Plasmas 21 (5) (2014) 055907.[41] J. Chen, X. Xu, Y. Lei, Extension of landau-fluid closure to weakly collisional plasma regime, Computer Physics Communications 236 (2019)128–134.[42] X. Xu, N. Li, Z. Li, B. Chen, T. Xia, T. Tang, B. Zhu, V. Chan, Simulations of tokamak boundary plasma turbulence transport in setting thedivertor heat flux width, Nuclear Fusion 59 (12) (2019) 126039.[43] F. Halpern, P. Ricci, S. Jolliet, J. Loizu, J. Morales, A. Mosetto, F. Musil, F. Riva, T.-M. Tran, C. Wersal, The gbs code for tokamak scrape-o ff layer simulations, Journal of Computational Physics 315 (2016) 388–408.[44] M. Francisquez, B. Zhu, B. N. Rogers, Multigrid treatment of implicit continuum di ff usion, Computer Physics Communications 236 (2019)104–117.[45] W. Zholobenko, T. A. Body, P. Manz, A. Stegmeir, B. Zhu, M. Griener, G. D. Conway, D. P. Coster, F. Jenko, Electric field and turbulence inglobal braginskii simulations across the asdex upgrade edge and scrape-o ff layer, Plasma Physics and Controlled Fusion (2021).[46] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, D. Karpeyev,D. Kaushik, M. G. Knepley, D. A. May, L. C. McInnes, R. T. Mills, T. Munson, K. Rupp, P. Sanan, B. F. Smith, S. Zampini, H. Zhang, PETScWeb page.URL [47] R. Falgout, HYPRE: High Performance Preconditioners Web page, Lawrence Livermore National Laboratory.URL [48] B. D. Dudson, J. Leddy, Hermes: global plasma edge fluid turbulence simulations, Plasma Physics and Controlled Fusion 59 (5) (2017)054010.[49] N. Li, X. Xu, T. D. Rognlien, B. Gui, J. Sun, D. Wang, Calculation of two-dimension radial electric field in boundary plasmas by usingbout ++ , Computer Physics Communications 228 (2018) 69–82.[50] C. Ma, B. Zhu, X.-Q. Xu, W. Wang, Machine learning surrogate models for landau fluid closure, Physics of Plasmas 27 (4) (2020) 042502., Computer Physics Communications 228 (2018) 69–82.[50] C. Ma, B. Zhu, X.-Q. Xu, W. Wang, Machine learning surrogate models for landau fluid closure, Physics of Plasmas 27 (4) (2020) 042502.