{'fraction_non_alphanumeric': 0.027949542418995794, 'fraction_numerical': 0.005243630967103636, 'mean_word_length': 4.604657610202384, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A tight-coupling scheme sharing minimum information across a spatial interface between gyrokinetic turbulence codes

A new scheme that tightly couples kinetic turbulence codes across a spatial interface is introduced. This scheme evolves from considerations of competing strategies and down-selection. It is found that the use of a composite kinetic distribution function and fields with global boundary conditions as if the coupled code were one, makes the coupling problem tractable. In contrast, coupling the two solutions from each code across the overlap region is found to be more difficult due to numerical dephasing of the turbulent solutions between two solvers. Another advantage of the new scheme is that the data movement can be limited to the 3D fluid quantities, instead of higher dimensional kinetic information, which is computationally more efficient for large scale simulations on leadership class computers.


I. INTRODUCTION
Fusion research will advance to the burning plasma stage with the planned completion and operation of ITER in the next decade. Guided by the understanding obtained from several fusion experiments as well as theory and simulation activities in the world, ITER is expected to attain ten-fold energy gain and realize plasmas that are well beyond the plasma regimes accessible in present and past experiments. It is thus critical for theory and simulation to be reliably predictive so that ITER can achieve its full potential through well-confined and disruption-free operation scenarios. Since fusion plasmas are complex entities, it is not sufficient to understand them piecemeal. The primary focus of this paper is on the development of the first steps that must be taken if we are to realize a Whole Device Model (WDM) of a fusion plasma, based on first principles, that will enable us to understand and predict the performance of existing fusion devices and next step facilities such as ITER.
In 2016, the Department of Energy (DOE) Exascale Computing Program (ECP) launched an ambitious simulation project called "High-Fidelity Whole Device Modeling of Magnetically Confined Fusion Plasma", simulating the whole plasma based on first-principles physics from the magnetic axis to the material wall. This is an ambitious challenge, because such a fusion plasma is characterized by multiple spatial and time scales. As a first step, one turns to gyrokinetic theory as a suitable point of departure, except that even in that case, the numerical methods that have been found to be most effective for the core and the edge plasma are different. This is because the underlying physics and plasma properties in these two domains of a turbulent tokamak plasma are quite different. Instead of attempting to describe the entire plasma by a single kinetic code, we have taken the approach of coupling two kinetic codes, one that works optimally in the core plasma, and the other in the edge (or boundary) plasma.
The first phase goal in ECP-WDM has been defined to be a tight kinetic coupling of core and edge micro-turbulence codes. The core and edge turbulence interact and affect each other, leading to possible emergent behavior that cannot be modeled accurately using one or the other code specific. The coupling between the two codes must be tight in order to yield an accurate and self-consistent global solution.
As a first step, we begin here by reporting the results from tests of the spatial coupling technique obtained by splitting the global gyrokinetic code XGC [1] into two executables: one for the core region and the other for the edge region. The strategy of using the same code for two executables located at different spatial regions, and treating them as separate codes, permits us to study the coupling techniques without the possible complication from any differences in the form of the governing equations used, algorithms used, and error from grid interpolation. Communication between the two executables is carried out using the ADIOS software [2].
In the present study, the plasma is modeled with the electrostatic gyrokinetic equations [3] in realistic tokamak geometry. In the gyrokinetic equations, particle motion in 6-dimensional (6D) phase space is reduced to 5D motion by analytically averaging out the gyrophase which represents fast gyration in a strong magnetic field. The target problem is further simplified by using the so-called "delta-f" scheme, in an effort to identify the most basic issues in the turbulence coupling through a spatial interface. In the "delta-f" scheme, the particle distribution function is split into two components f = f 0 + δf , with f 0 a Maxwellian background, consistent with an experimental plasma density, temperature and mean flow profiles, and δf a perturbation from f 0 in the particle distribution function responsible for micro-turbulence. In this study, the δf perturbation is solved with the Particle-In-Cell (PIC) technique. Collision and source terms are not included in the present simplified study. The plasma species modeled here are gyrokinetic ions and adiabatic (Boltzmann) electrons.

II. COUPLED EQUATIONS
In the system of coupled delta-f core and edge codes, a common background f 0 is employed in accordance with the specified radial density and temperature profiles. The perturbed particle distribution functions in the core and edge, respectively δf Core and δf Edge , are combined into a composite one where ̟ = ̟(ψ gc ) is a continuous connection function varying between 0 and 1 with respect to the gyrocenter radial position ψ gc . The radial domain is split into three regions: core (̟ = 1), overlap (0 < ̟ < 1), and edge (̟ = 0); see Fig. 1 and the description in the next sub-section. We use a simple linear form of ̟ here. The effects of a higher-order connection function will be studied later. The electrostatic potentialφ is obtained by substituting the composite distribution function δf , in the global electrostatic Poisson equation for the whole plasma volume as if the coupled system were a single code, The gyrokinetic Poisson operator, L, has coefficients that are dependent on the moments of f 0 , i.e., plasma density and temperature; see detailed equations in reference [4,5]. The source on the right-hand-side is a combination of core and edge charge density contributions that is calculated from the composite δf , where X is the gyrocenter position vector, ρ is the Larmor vector, x is the particle position vector in the configuration space, v is the particle speed scalar parallel to the magnetic field line, µ is the magnetic moment, α is the gyroangle, and δ[· · ·] is the Dirac operator. For consistency, it is crucial to have ̟ defined in gyrocenter space, ̟ = ̟(X), before applying the guiding-center pull-back operation, δ[X + ρ − x]. An important consequence of Eq. (3) is that the kinetic coupling is enabled by kinetic integration of the code-local distribution functions δf C and δf E in each code, without the necessity for communicating the 5D information between the codes. This has a significant advantage in large-scale computing where the data movement is expensive, because only the 3D field datan C andn E need to be communicated between the codes for solving the field equation consistently with kinetic physics. The unified, self-consistent fieldφ is then used, in both the core and edge codes, to evolve the distribution function on each side in accordance with the gyrokinetic Vlasov equation Here q the electric charge, m the particle mass, α the gyroaveraging operation, b = B/B with B the magnetic field and B its magnitude, and B ⋆ the effective magnetic field; see Ref. [4,5]. In this core-edge coupling scheme, the only boundary condition to be imposed on the coupled system is the particle and solver boundary condition at the outer boundary of the edge code, as if the coupled system is from a single whole-volume code. For example, the plasma particles are absorbed at the material wall and recycled back into the plasma region as neutral particles, and the electrostatic potential is grounded at the material wall. Importance of enforcing this self-consistency between the unified global fieldφ and the particle distribution function in obtaining the correct turbulence solution across the interface region cannot be overstated. As is to be shown later, coupling turbulence solutions through an interface using solutions from each code turns out to be problematic.
In this scheme, the time-stepping of the global core and edge distribution functions is achieved by pushing the composite distribution function independently in each code, but using the common global potential field solution. Injecting δf as defined in Eq. (1) into the gyrokinetic total derivative operation leads to If δf C and δf E are the same at t = 0 and are evolved using the same gyrokinetic equation, including the coupled self-consistent fieldφ, then δf C ≃ δf E locally in the interface region at all times t > 0, and Numerical errors may not guarantee that δf C = δf E is exactly satisfied, in particular for PIC codes where there might be a systematic sampling difference between δf C and δf E . Enlarging the overlap region decreases ∇̟ and increasing the number of markers decreases the sampling difference on δf C − δf E . This error source in the coupling scheme needs to be carefully examined through convergence study. The tight coupling scheme based on these equations is applied at each stage of XGC's Runge-Kutta time integrator with a four-step algorithm: 1. Prior to solving for the field, both the edge and core simulations compute their contribution to the charge density, Eq. (3).
2. The core simulation sends its charge density contribution to the edge simulation.
3. The edge solves the gyrokinetic Poisson equation globally including the core region, Eq. (2), and sends the fielď φ to the core code. We use the edge solver here because, the edge code is usually more general and can cover the whole volume. The field solver could be a separate routine. Using only one global field solver guarantees that the exact same field is used by both codes.
4. Both simulations push their δf with the same gyrokinetic Vlasov equation, by using exactly the same global fieldφ distribution, Eq. (4).

III. RADIAL DOMAIN DECOMPOSITION WITH OVERLAP AND BUFFER REGIONS
The interface between the core and edge codes is an overlap region where 0 < ̟ < 1. As detailed in the previous section, both codes consistently contribute in this overlap region to the Poisson equation's source term defined by Eq. (3) and consistently respond to the Poisson-solved field through Eq. (4).
In addition to the overlap region, each of the core and the edge simulation modules have buffer zones. The buffer zones for a particular module is spatially contiguous with the overlap region, and extends into the spatial region for the other module. The buffer region for a particular module accommodates the particles that travel outside of the interface region into the other region, and return. The existence of these buffer regions also ensure continuity of the turbulence fluxes at the inner and outer edges of the overlap region. In these regions, marker particles are evolved with the coupled Eq. (4) but do not contribute to the Poisson equation source term, Eq. (3). In the case shown in Fig. 1, the buffer region of the core module is located between the outer edge of the interface region (ψ = 0.76) and 0.95, where ̟ = 0. In the present study, we have chosen the buffer region of the edge simulation to extend all the way to the magnetic axis from the inner edge of the overlap region, but it can be chosen to have a more limited spatial extent. Since the plasma volume decreases with r 2 in a torus, with r being the minor radius, extension of the inner buffer region all the way to the magnetic axis does not increase significantly the computational cost but increases the accuracy of the coupling. For the purpose of the present study, including only the outer buffer region will be sufficient for the study of the effect of this region on code coupling.
The outer buffer region for the core code cannot be extended all the way to the material wall because the equations underlying the core code become invalid near the magnetic separatrix surface and in the scrape-off layer. Thus, a finite radial width is needed, separating the outer buffer region from the magnetic separatrix. If a particle from the core region exits the outer buffer region or equivalently, exits the simulation region of the core code into the edge region, it is re-injected back to the buffer region along the gyrocenter orbital path in a way which preserves density, momentum, and energy for each marker particle.

IV. COUPLED SIMULATION WITH A FIELD SOLVER
We demonstrate the fidelity of our coupling scheme, by comparing in Fig. 2 the results obtained with a reference simulation (a), a coupled simulation (b) and the time-integrated heat-flux across the overlap region in Fig. 3. Specifically, in these figures we choose to use the ion heat flux, with S defines surface averaging. Heat flux is one of the main quantities of great interest in measuring the performance of a magnetically confined fusion plasma. The geometry of our equilibrium is a DIII-D-like plasma including an X-point. The ion temperature profile is chosen in such a way that it produces turbulence in the edge region, with a strong temperature gradient just outside of the overlap region. It is well-known that turbulence then propagates radially inward [6]. There is no heat source in the simulation. Thus, the temperature gradient relaxes as a result of the turbulent transport, and the turbulence eventually decays. For the boundary condition of the core-edge coupled simulation, the electrostatic potential is zero and particles are absorbed at the material wall outside the magnetic separatrix surface. For the present simulation, the material wall is placed along a surface that conforms to a magnetic surface at ψ n = 1.11.
In Fig. 2 subplot (b), it is shown that the turbulence in the coupled simulation can propagate freely from the edge to the core domain, i.e., from one executable to the other, with a high level of fidelity for the problem. To be more quantitative, as measured in Fig. 3 at the end of the simulation, there is only a 3% difference in the time-accumulated turbulent heat flux, dtQ i (ψ; t), passing through a surface in the overlap region, which is the most sensitive region for the coupling challenge. For perspective, it is worth keeping in mind that experimental error bars, against which we eventually need to validate the simulation results, are usually greater than 10-20%.
We remark on a notable result of the present study that may appear counter-intuitive. It may appear upon first glance that the most straightforward mathematical way to couple two codes is by having each code use its own solvers and outer boundary conditions, with the solutions matched at the inner interface. Our finding is that this is not so, and that the spatial coupling becomes much easier if two codes solve for the same composite solution, subject to the boundary conditions as if the two codes were one. As an example, we present here a case which uses separate solvers, with the boundary condition for the core code at the outer end of the buffer region ψ n = 0.95. Since ψ n = 0.95 is located multiple radial-orbit-widths and turbulence-correlation-lengths away from the core region ψ n < 0.76, the imposition of a separate solver boundary condition for the core code does not affect significantly the coupled solution. The result of this case is shown in the subplots (c) of Fig. 2, and the red line in Fig. 3, which is from the following two boundary conditions for the solvers executed in each codes: As can be seen from Fig. 2 (c), turbulence is partially suppressed in the overlap region, and the heat flux through it is much lower than the reference case (Fig. 3). A detailed study reveals that this is caused by the phase mismatch between the core and edge solutions of the fluctuating turbulent potential field φ C and φ E , resulting in a partial cancellation between them. An artificial solver boundary condition at ψ n = 0.95 is causing this phase mismatch, by not allowing for the phase information to be correctly communicated between the core and edge regions.

V. INVESTIGATION ON THE INFLUENCE OF THE OVERLAP WIDTH IN A SIMPLIFIED GEOMETRY
Another significant finding from our study is that the overlap width, w, must be larger not only than the radial orbital width but also the radial turbulence correlation length, λ r . To study this problem, we carry out a nonlinear simulation scan with respect to the overlap width, from w = 0 × λ r to w ≃ 3.4 × λ r . λ r is evaluated to be 2.3 cm in the overlap region from the reference simulation. To avoid possible complications from the flux-surface shape, we use a circular cross-section torus, called "cyclone-base-case geometry" with the temperature gradient R/L Ti = 6.9 and the density gradient R/L N = 2.2. No heat source is applied, such that the plasma is left free to relax through the growth of δf via the turbulent heat flux down the gradient. For this relaxation problem, we measure the ion temperature gradient R/L Ti , plotted in the first row of Fig. 4, as well as the radial average of the ion heat flux and the ion temperature gradient in the second row. The ion temperature gradient relaxation is the result of the ion heat flux integrated over time and is a very good indicator of the accuracy of the coupling. Towards the end of the simulation, since the temperature gradient has relaxed and the turbulence is weak, the comparison is less meaningful.
The reference and the coupled simulations agree well when the overlap region width, w, is wide enough. In our scan, the best agreement is obtained when w is a few times the local correlation length of turbulence, λ r . Even when w ≃ 1.7λ r , there is a good agreement between the coupled (green) and the reference (red) cases in the radial-averaged heat flux and temperature gradient in subplots (f) and (g), respectively,. However, the zonal structure on R/L Ti near r/a ≃ 0.7 (red arrows) of the reference simulation (a) shows better agreement with the results with (e) w ≃ 3.4λ r than with (d) w ≃ 1.7λ r .

VI. CONCLUSION
A new tight-coupling scheme is introduced that enables spatial coupling of large-scale kinetic turbulence codes across an interface. Several key points of this scheme are explained. Distribution functions from the core and edge codes are combined into one by using a connection function at the overlap region. Using this connection function, the charge density is calculated from the kinetic information obtained from each code and communicated to the solver routine, whereby the two 3D charge density distributions are combined into one to determine the source term on the right hand side of the global Poisson equation. This method allows us to avoid expensive communication of large-scale 5-dimensional kinetic data. For successful coupling, the key is that the global Poisson equation and the boundary condition needed to solve it are as if the coupled codes are one. Furthermore, it is found that the overlap region at the coupling interface needs to be of similar size or wider than the radial correlation length of turbulence and radial orbit width.
In related work, using simplified model equations, Ricketson et al. [7] have recently shown that classical additive Schwarz-type methods hold considerable promise for addressing coupling problems of the kind considered in this paper. They obtain analytical as well as numerical results on convergence for their model problem. It would be worthwhile to consider using similar techniques why the method proposed in this paper is effective.
Spatial code coupling is an active research area. As the next step in the ECP WDM project, the present scheme is being successfully extended to coupling two different core-edge gyrokinetic codes (GENE [8] in the core and XGC in the edge), which involves grid mapping as an extra ingredient. This work will be reported in the near future.