第 21 回
核燃焼プラズマ
統合コード研究会
21st Burning Plasma Simulation
Initiative (BPSI) Meeting
日時：2023 年 12 月 21 日(木)－22 日(金)場所：九州大学筑紫キャンパス応用力学研究所
2 階大会議室 および オンライン(Ver. 1, 2024 年 3 月 1 日発行)

合同会合第 21 回核燃焼プラズマ統合コード研究会
21st Burning Plasma Simulation Initiative (BPSI) Meetingおよび
核融合エネルギーフォーラムサブクラスター会合
シミュレーションクラスター 炉心プラズマモデリングサブクラスター
プラズマ物理クラスター 定常運転・制御サブクラスター(Ver.2)
日時：2023 年 12 月 21 日(木)－22 日(金)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室 および オンライン

12 月 21 日(木)(特定研究「計測・シミュレーション・モデリングを組み合わせた統合診断」研究集会との合同会合)(座長：糟谷)9:45 – 10:00 講演 D-1
糟谷（九大）
特定研究「計測・シミュレーション・モデリングを組み合わせた統合診断」の実施
状況
10:00 – 10:30 講演 D-2
古賀（兵庫県立大）Research on a reconstruction method of micro-wave imaging using machinelearning
機械学習を用いたマイクロ波イメージング再構成手法の研究
10:30 – 11:00 講演 D-3
佐藤(雅)（核融合研）
Simulation data analysis of nonlinear structures of MHD instabilities in toroidal
plasmasトロイダルプラズマにおける MHD 不安定性の非線形構造のシミュレーションデー
タ解析
11:00 – 11:30 講演 D-4
沼波（核融合研）
Development of a platform for quantitative understanding of plasma profile
formation mechanisms using first-principal simulation, data science, andnumerical measurement
プラズマ分布形成の定量的理解に向けた第一原理シミュレーションとデータ科学、
数値計測による基盤開発
11:30 – 12:00 糟谷（九大）
discussion 議論（一部関係者のみ）
12:00 – 13:00 昼休み
113:00 – 13:20 会議登録 Registration
13:20 – 13:30 事務連絡 Business announcement(座長：村上)13:30 – 13:55 講演 1-1
本多 (京大)Assessment and feasibility study of core plasma performance for JA-DEMO
using GOTRESS+
13:55 – 14:20 講演 1-2
福山 (京大)Progress of integrated code TASK and the development of interface with IMAS
14:20 – 14:40 講演 1-3
持永 (九大)Effect of toroidal rotation on impurity transport in tokamak improved
confinement
14:40 – 14:55 講演 1-4
田原 (京大)The effect of beam-beam coulomb collisions on the energetic particle distribution
in LHD
14:55 – 15:15 集合写真 group photo、休憩 break(座長：福山)15:15 – 15:40 講演 1-5
矢木 (量研)GPU acceleration of RMHD model and application to particle transport inplasma peripheral
15:40 – 16:05 講演 1-6
奴賀 (核融合研) online
Degradation of fast-ion confinement depending on the neutral beam power
without EP-driven instability in LHD
16:05 – 16:30 講演 1-7
森下 (京大)Development of an adaptive predictive control system using data assimilation
at LHD
16:30 – 16:45 講演 1-8
佐藤(真) (総研大)Study of model based control for tokamak plasma profile
16:45 – 17:00 講演 1-9
鈴木 (総研大)
Prediction of radiative collapse by analyzing videos of plasma discharge with
CNN in LHD
17:00 散会2

12 月 22 日(金)9:00 – 10:30 ポスター講演 Poster(online, 開室 start 8:50、閉室 end 12:50)
10:30 – 10:40 休憩(座長：糟谷)10:40 – 11:05 講演 2-1
相羽 (量研) online
Identification of plasma conditions affecting MHD stability in QH-mode andELMy H-mode plasmas in DIII-D
11:05 – 11:30 講演 2-2
佐々木 (日大) online
Interaction between turbulence and zonal flows in the presence of mean shearflow
11:30 – 11:45 講演 2-3
竹本 (名古屋大) online
Comparison of W spatial distribution and time evolution by impurity transport
calculation with soft X-ray measurement data
(Session Leader：村上)
11:45 – 12:05 議論 Discussion
12:05 – 12:10 事務連絡 Business announcement
12:10 – 13:30 昼休み Lunch break(ワークショップ「不純物輸送」 Lunch time workshop “Impurity transport”
12:15 – 13:00 林 (QST)、藤田 (名大)、糟谷 (九大) online
昨年に続く共同研究作業会を昼休みに開催します、関係者のみによる集まりです)(核融合エネルギーフォーラムサブクラスターとの合同会合)
13:30 – 13:40 林（QST）
連絡事項
13:40 – 14:10 宮戸（QST）
IFERC 計算機シミュレーションセンターの現状報告
14:10 – 14:40 横山（核融合研）
ITER に於ける統合モデリング活動の報告
14:40 – 15:10 若月（QST）
ITPA 統合運転シナリオグループ活動報告
15:10 – 15:30 星野（慶應大）炉心プラズマ＋炉工学炉材料モデリングサブクラスター合同会合の報告
15:30 – 16:00 林（QST）
今後の予定の議論
16:00 散会
3

Poster number
21th BPSI meeting (2023.12)(Ver.1)(コアタイム 前半: 奇数番、後半: 偶数番)
PA-1 村上定義 (京大) S. Murakami
Modeling of ECH supra-thermal electron transport in the presence of magnetic andelectrostatic fluctuations
PA-2 成田絵美 (京大) E. Narita
Introduction of hydrogen isotope effects to a turbulent transport model DeKANIS and
its application to integrated simulations
PA-3 登田慎一郎 (核融合研) S. Toda
Research of turbulent transport due to dissipative trapped electron mode in tokamak
plasmasPA-4 糟谷直宏 (九大) N. Kasuya
Development of integrated turbulence diagnostic simulator and its application to torus
device measurement
PA-5 福永龍平 (京大) R. Fukunaga
Verification of charge separation model in generation of radial current with multi-fluid
transport code TASK/TX
PA-6 古田原拓実 (日大) T. Kodahara
Extraction of abrupt turbulence driven particle transport using multi-field SVD
PA-7 可児和寿 (名大) K. Kani
Optimizations of NB deposition calculations in the integrated code
PA-8 桑宮幸佑 (九大) K. Kuwamiya
Simulation of PCI experimental diagnostic for density fluctuations in tokamak
plasmasPA-9 宮本琉耶 (九大) R. Miyamoto
Evaluation of nonlinear mechanism to induce inward turbulent particle fluxes intokamaks
PA-10 小山一輝 (九大) K. Oyama
Analysis of parallel flow shear driven mode including effect of finite ion temperature
gradient4「計測・シミュレーション・モデリングを組み合わせた統合診断」
特定研究の実施状況
Present status of designated joint research “Synthetic diagnostics
combining measurement, simulation and modeling”
N. Kasuya
RIAM, Kyushu Univ.

九州大学応用力学研究所では「計測・シミュレーション・モデリングを組み合わせた統合診断」をテーマとした特定共同研究を実施している。本特定研究では、様々な分
野で行われている実験計測、複合シミュレーション、データ駆動モデリングを組み合わせた統合的診断手法の共通性や応用性を探ることで、液体・気体・プラズマのミク
ロ・マクロ的様相を理解するためのさらなる手法の進展を目指している。
以下の 6 件のサブテーマを実施している。
・データ駆動科学的手法を用いた非平衡開放系の時空間ダイナミクス研究 (研究代表者 日本大学 佐々木真准教授)
・計測との比較に向けた不純物輸送の統合モデリング (研究代表者 京都大学 本多充教授)
・機械学習を用いたマイクロ波イメージング再構成手法の研究 (研究代表者 兵庫県立大学 古賀麻由子准教授)
・トロイダルプラズマにおける MHD 不安定性の非線形構造のシミュレーションデー
タ解析 (研究代表者 核融合科学研究所 佐藤雅彦助教)
・直線磁化プラズマにおける乱流構造の解析 (研究代表者 九州大学 山田琢磨教授)
・プラズマ分布形成の定量的理解に向けた第一原理シミュレーションとデータ科学、
数値計測による基盤開発 (研究代表者 核融合科学研究所 沼波政倫教授)
本セッションではそのうちより 3 件の研究進捗状況を発表いただき、相互理解を深めるとともに、共同研究進展の可能性について議論する。

Research on a reconstruction method of micro-wave imaging
using machine learning
Mayuko Koga
Graduate School of Engineering, University of Hyogo, Himeji, Hyogo 671-2280, Japan
Plasma imaging plays a pivotal role in fusion plasma study. Sub-microwave and microwave
diagnostic techniques, such as reflectometry, offer the advantage of straightforwardidentification of measured points [1]. To overcome the challenge of achieving high resolution
imaging with a limited number of antennas, there is a need to explore method for image
reconstruction from sparse measurements. Tsuchiya et al. previously demonstrated thefeasibility of reconstructing simple images from complex-amplitude microwave data using
Tikhonov-Phillips regularization techniques [2]. However, generating videos from numerous
images remains impractical due to the extensive computational time required for each image
reconstruction. In this paper, perform the principal verification of reconstructing object surface
images by using deep learning.
In various fields, the significance of the quality and scale of training data is well-recognized
for supervised learning. Ideally, training data should be sourced from actual plasma experiments.
However, due to the inherent difficulty of controlling the cutoff surface pattern in plasma
experiments, we undertook numerical calculations of reflected microwave distribution across awide range of patterns. Subsequently, this dataset, consisting of pattern-distribution pairs, wasemployed for our training purposes. We made objects with a protrusion. The footprint of the
protrusion is set to elliptical, whose length of major and minor axis are randomly set in therange of 1 mm to 65 mm, respectively. The protrusion shape is a Gaussian function like, whose
top is 1 mm height. The position and the rotation angle of the protrusion are also randomly set.
The reflectance of the object is set to 1. Thirty gigahertz microwave is irradiated to the object
and completely reflected at the surface of the object. The reflected microwave signals arereceived by an antenna array. The antenna array consists of horn antennas arranged in 16 × 16
pieces, with each antenna having a size of 10 mm square. We randomly changed the structureof protrusion and calculated the corresponding received microwave signals.
We used convolutional neural network (CNN) as the neural network for deep learning. We
tried three patterns as input channels, absolute values of complex amplitudes, absolute valuesof complex amplitudes and phase, real and imaginary part of complex amplitudes. The imageof received microwave signals is input into the network with 16 × 16 pixels. It passes througha convolutional layer (karnel size: 6 × 6, number of filters: 50) and a max pooling layer (window

size: 2 × 2). After that, it passes through 3 layers of the fully connected layer which contains
50 perceptron and one reshape layer. Tanh function was used as the activation function. After
passing through the reshape layer, a predictive image of object surface is outputted with 128
× 128 pixels. This means that the output image is upconverted from the input image.

Figure 2 and 3 show the example of reconstruction. The input channels take absolute valuesof complex amplitudes (Fig. 2) and real part and imaginary part of complex amplitudes (Fig.
3). The color bars mean normalized intensities in (a), surface positions in the depth directionin (b) and (c). The numbers above images show the numbers of training data used for eachmodel. Both figures show that reconstructed images capture features of object images. As thenumber of training data increases, it is clear that small structures of objects are well
reproduced. It was found that two types of input data (real-imaginary part) are better than
one type of input data (absolute values). Reconstruction calculations on a trained model takes
about milliseconds computing time, and could be applied for real-time feedback, etc in future.

References[1] Y. Nagayama, N. Ito, D. Kuwahara, H. Tsuchiya, and S. Yamaguchi, Rev. Sci. Instrum. 88,
044703 (2017).
[2] H. Tsuchiya, N. Iwama, S. Yamaguchi, R. Takenaka and M. Koga, Plasma Fusion Res. 14,3402146 (2019).

Detected Object Predicted image
signal image （Number shows the number of training data）

Fig. 1. (a) Received microwave signals, (b) object surface images and (c) corresponding
reconstructed surface images by models with different numbers of training data. The
input channels take absolute values of complex amplitudes. The numbers above
images show the numbers of training data.
Detected Object Predicted image
signal image （Number shows the number of training data）

Fig. 2. (a) Received microwave signals, (b) object surface images and (c) corresponding
reconstructed surface images by models with different numbers of training data. The
input channels take real part and imaginary part of complex amplitudes. The
numbers above images show the numbers of training data.

Simulation data analysis of nonlinear structures of MHD instabilities in toroidal plasmas(Analysis of infernal modes in tokamak plasmas)
M. Sato1, Y. Todo1, N. Aiba2 and M. Takechi2
1
National Institute for Fusion Science, Toki, Gifu, 509-5929, Japan
2
National Institutes for Quantum Science and Technology, Naka, Ibaraki, 311-0193, Japan
In order to realize a steady-state fusion reactor, high beta plasma must be stably confined fora long time. In high beta plasmas, the bootstrap current, which is a non-inductive current,
increases and discharges with an increased proportion of bootstrap current are one of the steadystate operation scenarios in tokamak plasmas. Since the bootstrap current takes its maximum
value away from the magnetic axis, resulting in a flat q-profile ornegative magnetic shear configuration, the magnetic shear in thecore region becomes weak. It has been pointed out that the presenceof a large pressure gradient in such a weak magnetic shear region
can lead to infernal modes with low toroidal mode number[1]. In[2, 3], the linear ideal MHD analyses showed that the occurrenceof poloidal beta collapse in JT-60U is consistent with the violationof the stability boundary of the ideal infernal modes. However,since the previous studies are based on the linearized ideal MHD
model, the influence of the saturated state of both the ideal infernal Fig.1 Profiles of pressure (solidmodes and the resistive infernal modes on the plasma confinement curve) and the safety factor
has not been clarified yet. For the highly accurate evaluation of the (dashed curve) for the MHD
equilibria analyzed in this study.
beta limit, not only linear stability but also nonlinear saturated stateshould be investigated for both the ideal and the resistive modes.
For the evaluation of the beta limit, the effects of kinetic thermal
ions might be important since the linear growth rate of theinstabilities becomes small for the plasmas near the beta limit. For
example, in the LHD plasmas, it was found that the effects of the
kinetic thermal ions play an important role for the stability beta
limit[4]. In this study, kinetic-MHD hybrid simulation analysis forthe infernal modes is carried out in order to clarify the effects of the
kinetic thermal ions on the linear stability and the saturated statesof the infernal modes.
Fig.2 Dependence of the linear
The numerical calculations are performed by MIPS code growth of the n=3 mode on the(MHD model) and MEGA code (Kinetic-MHD hybrid model) In magnetic Reynolds number
the MEGA code, thermal ions are treated by the drift kinetic model (S) for β0 = 5.7%.
and the electrons are treated by the fluid model[4, 5]. The periodic
boundary condition is imposed at φ= 0 and φ = 2π/3 so that
only the n = 3 mode and its harmonics are included where φis thetoroidal angle. In this study, two MHD equilibria with circularpoloidal cross section are analyzed. Figure 1 shows the profiles ofthe pressure and the safety factor of the two equilibria. The central
beta value β0 of the two equilibria is β0 = 5.7% and β0 = 4.1%. For
β0 =5.7%, an n=3 ideal infernal mode is unstable where n is thetoroidal mode number. On the other hand, for β0 = 4.1%, the n=3
ideal infernal mode is stable while the n=3 resistive infernal mode
is unstable.
Fig.3 Dependence of the linear
Figure 2 shows the dependence of the linear growth rate on growth of the n=3 mode on thethe magnetic Reynolds number (S) for β0 = 5.7%. Although the magnetic Reynolds number(S) for β0 = 4.1%.

Fig.4 Time evolution of the pressure profile on apoloidal cross section obtained from the MHD
model (the upper rows) and the kinetic-MHD
hybrid model (the lower rows) for (β0, S) =(5.7%, 108).

Fig.5 Time evolution of the pressure profile on apoloidal cross section obtained from the MHD
model (the upper rows) and the kinetic-MHD
hybrid model (the lower rows) for (β0, S) =(4.1%, 107).

linear growth rate depends on the S number in low S regime, the linear growth rate becomes
independent on the S number in high S regime, which indicates that the unstable mode is an
ideal mode. The kinetic-MHD hybrid simulations with kinetic thermal ions were conducted for
S=105 and S=108, with their linear growth rates illustrated by the red circles. For S=105, thelinear growth rate obtained from the kinetic-MHD hybrid model with kinetic thermal ions isalmost the same as that obtained from the MHD model without kinetic thermal ions. When
S=108, the linear growth rate obtained from the kinetic-MHD hybrid model with kinetic thermal
ions is smaller than that obtained from the MHD model without kinetic thermal ions. The linear
growth rates for β0 = 4.1% obtained from the MIPS code and the MEGA code are shown in Fig.
3. The linear growth rate obtained from the MHD model without kinetic thermal ions decreases
as the S number increases, indicating a resistive infernal mode remains unstable. For S=107, thelinear growth rate obtained from the kinetic-MHD hybrid model with kinetic thermal ions issmaller than that for the MHD model without kinetic thermal ions.
Figure 4 shows the time evolution of the pressure profile on a poloidal cross section for(β0, S) = (5.7%, 108). For both the models, at the saturated state, the influence of the instabilities
expands to the almost the whole region and the pressure profile is significantly flattened. Thus,
the influence of the effects of the kinetic thermal ions is weak for the ideal infernal modes.
Figure 5 shows the time evolution of pressure profile on a poloidal cross section for (β0, S) =(4.1%, 107). The pressure profile obtained from the MHD model is significantly flattened at thesaturated state. On the other hand, for the kinetic-MHD hybrid model, the pressure profile is
not flattened at the saturated state. This indicates that kinetic thermal ions play an essential role
for the suppression of pressure profile flattening due to slowly growing resistive MHD
instabilities.
[1] J. Manickam et al, Nucl. Fusion 27 (1987) 1461.
[2] T. Ozeki et al., Nucl. Fusion 33 (1993) 1025.
[3] T. Ozeki et al., Nucl. Fusion 35 (1995) 861.
[4] M. Sato and Y. Todo, Nucl. Fusion 61 (2021) 116012.
[5] Y. Todo et al., Plasma Phys. Control. Fusion 63 (2021) 075018.

Assessment and feasibility study of core plasma performance
for JA-DEMO using GOTRESS+
M. Honda, N. Aiba1, M. Furukawa2, Y. Masamoto2, and R. Otani
Kyoto University, Nishikyo-ku, Kyoto 615-8530 Japan
1

National Institutes for Quantum Science and Technology, Rokkasho, Aomori 039-3212 Japan2

Tottori University, Tottori, Tottori 680-8550 JapanGOTRESS+ extension for burning plasmas
In DEMO, a plasma is alpha-heated almost exclusively in a stationary state. Given
that auxiliary heating systems have functioned as actuators to control the plasma, inburning plasmas, where the contribution of auxiliary heating to the total heating power isquite small, it is not always obvious whether and how a stationary state can be achieved
and then be stably maintained. Therefore, it is important to verify them through numericalsimulations. Recent updates to convergence algorithms in GOTRESS enable the robust
identification of stationary states while accounting for alpha heating and large radiation
loss channels, both of which are unique to burning plasmas. GOTRESS+, a novel
integrated model GOTRESS+ [1] as illustrated in Fig. 1, is therefore ready to performpredictive burning plasma simulations. We verify whether the target values for the JADEMO plasma [2] can be achieved with GOTRESS+. The plasma confinement
performance is greatly affected by the choice of turbulent transport model. In this case,
the CDBM turbulent transport model [3] was selected.

Assessment of JA-DEMO operation scenario with ECH only
Electron cyclotron heating (ECH) can realize a localized absorption profile,allowing adjustments to the current profile shape to reduce magnetic shear locally andpromote internal transport barrier (ITB) formation. The research investigates whether it
is possible to achieve the target values with less power by adjusting the injection
parameters of ECRF waves. For simplicity, it assumes that the ECH power is injectedfrom a single launcher.

Figure 2 shows the summary of the GOTRESS+ simulation results for the scenario
with ECH only. In this case, O-mode wave with 185GHz and 90MW was injected fromthe mirror with toroidal angle of 150° and poloidal angle of 25° located at(R,Z)=(11.5,0.9) in meters, as shown in Fig. 2 (a). As evident from Fig. 2 (d), the localized
absorption of the EC wave results in 2.9 MA EC current generation. Consequently, as
observed in Fig. 2 (b) near the absorption peak, a negative magnetic shear is formed,
suppressing turbulent transport there and leading to the formation of the electron and ion
ITBs, as shown in Figs. 2 (f) and (g). In this case, all of the figures of merit meet the target
values [2]. However, we note that the global ideal MHD stability has not been achieved,
even if ion diamagnetic drift stabilization is taken into account. To achieve the target valueof Q=13, it is necessary to obtain a fusion output of at least 1.17GW with ECH power of
90MW, which inevitably requires high plasma confinement through ITB formation. This,in turn, necessitates a high β value, making MHD modes more prone to be unstable.
Among the target values, the requirement for the Q value seems the most stringent, andrelaxing this target value somewhat could allow for an MHD-stable scenario.

References[1] M. Honda et al., Nucl. Fusion 61, 116029 (2021).
[2] Y. Sakamoto et al., IAEA 27th Fusion Energy Conf., FIP/3-2, Gandhinagar (2018).
[3] M. Honda and A. Fukuyama, Nucl. Fusion 46, 580 (2006).

Fig. 1: The flowchart of a GOTRESS+ simulation. When not using EPED1, one proceeds the convergencecheck along the magenta dashed arrow extending from GOTRESS after finishing a GOTRESS simulation.
If the NB-driven current calculated by OFMC is used in the next ACCOME calculation as an option, thedriven current data is passed along the blue dashed arrow.

Fig. 2: GOTRESS+ prediction of stationary-state profiles for the JA-DEMO scenario with ECH only. (a)The equilibrium with EC ray and the resonant layer superimposed, (b) the safety factor q and themagnetic shear s profiles, and (c) density profiles are given, while profiles of (d) the heat source, the sink,(f) the temperature T and diffusivity c for electrons and (e)(g) those for ions are calculated

GPU acceleration of RMHD model and application to particle transport in plasma
peripheral
M. Yagi, H. Todoroki, N. Kasuya
National Institutes for Quantum and Radiological Science and Technology,
Rokkasho Fusion Institute, Rokkasho, Aomori, Japan
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka, Japan
1. Introduction
In fusion plasmas, the burn control is important issue to sustain the self-ignition condition.
To understand particle transport mechanism is essential to supply the fuel ions at plasma
peripheral. We are developing the transport-MHD model based on the reduced MHD (RMHD)
model to investigate this problem. In the present work, the implementation of R5F code (5-field
RMHD) on GPU is performed using OpenACC. It is found that the subroutine where nonlinear
coupling is calculated using the pseudo-spectral method is the most expensive part of the code.
We have reduced communications between CPU and GPU in between cuFFT calls. As the result,the elapse time of this subroutine is reduced more than half.
2. Numerical Scheme and Parallel Algorithms
The code is described on the cylindrical coordinates (𝑟, 𝜃, 𝜁), where the finite difference
method is applied for 𝑟 direction and the spectral method is applied for 𝜃 and 𝜁 directions.
For parallel algorithms, we have adopted the 1-dimensional domain decomposition method.
The 5-field RMHD model is schematically written as𝜕𝑈 = 𝐿(𝑈) + 𝑁(𝑈, 𝑈)
𝜕𝑡where 𝑈 is the vector and 𝐿, 𝑁 are the linear and nonlinear operators, respectively. The linear
terms are solved by Crank Nicolson implicit method and the nonlinear terms are evaluated by
the predictor-corrector method (namely, the first step is the explicit and the second step is theimplicit way). They are given by
𝑈 !"#/% − 𝑈 !
𝑈 !"#/% + 𝑈 !
= 𝐿2
3 + 𝑁(𝑈 ! , 𝑈 ! )∆𝑡/2
2
𝑈 !"# − 𝑈 !
𝑈 !"# + 𝑈 !
= 𝐿2
3 + 𝑁4𝑈 !"#/% , 𝑈 !"#/% 5∆𝑡2

The matrix inversion can be performed by the LAPACK library.

Figure 1: Main loop in R5F
The benchmark environment provides that Intel Xeon Gold 6248R (2.8GHz, 24cores) x2 for
CPU and NVIDIA V100 GPUs with Nvlink 2.0 and 32GBx4 memory. For the compiler,
NVIDIA HPC SDK 22.5 (nvfortran) and CUDA 11.4 is adopted. The compiler options are
given by -O4 -fast -mcmodel=medium -fastsse -I/fftwpath/include -Mpreprocess -DUSEMPI2
for CPU execution and -O4 -fast -mcmodel=medium -fastsse -acc -Mcuda I/cudadir/11.4/include -Mpreprocess -DCUFFT -DUSEMPI2

Figure 2: Cost evaluation for each

Figure 3: Comparison of execution

subroutine using CPU version

time for TMRHS between CPU andGPU version

Figure 1 shows the main loop in R5F and we found that the cost of subroutine TMRHS ismost expensive part shown in figure 2. As the first step, we simply replaced FFTW to cuFFT
which is executed on GPU. It becomes slower than FFTW on CPU shown in figure 3. It isfound that the elapse time for data transfer between CPU and GPU is dominant in FFT area for
GPU version. If we compare the elapse time of FFT execution, cuFFT is roughly 44 times fasterthan FFTW. Therefore, the next step is to suppress frequent data copy between CPU and GPU.
It can be done using enter data and exit data directive as well as kernels directive.

Figure 4: Enter data area is shown blue box.

The variables which are used by cuFFT can be stationed on GPU memory if we use enter data
directive such as!$acc enter data create(W1C(1:LYMED,1:KZMWD,S:E),…)!$acc host data use device(W1C(1:LYMWD,1:KZWD,S:E),…)DO I=S:E
istat=cufftExecZ2D(cuplan,W1C(1,1,I),W1R(1,1,I))
END DO

Figure 5: Tuning results the 4th bar (from left) shows the tuned resultfor 1GPU and the 5th bar the one for 4 GPUs.
The tuning results show that 1 GPU case is 1.81 times faster than CPU case and 4 GPU case is
3.91 times faster than CPU case. The key for GPU implementation is to suppress thecommunication between CPU and GPUs as much as possible. The off-load technique isessential for this purpose.
4. Summary and discussion
We have implemented 5-field RMHD model on GPU machine (SGI HPE8600). The goal isset that the elapse time of the subroutine TMRHS by GPU is more than double shorter over
CPU version. The communications between CPU and GPU before and after the function callof cuFFT library are reduced using enter data directive. As the result, we obtained fairly goodresult. Also, we extended the 4 GPU use in a node. In this case, GPU version gives almost 4
times faster than that of CPU version. The optimization of MPI communication using NVLink
for TMRHS should be performed. In addition, the subroutine GATHER should be implemented
on GPU. There are left for future work.
Acknowledgements
The computations were carried out on SGI HPE8600 supercomputer at QST and JAEA.
This work was supported in part by the Collaborative Research Program of ResearchInstitute for Applied Mechanics, Kyushu University (2023CR-FP-25).

Degradation of fast-ion confinement depending on the neutral beam power
without EP-driven instability in LHD
H. NUGA1 , R. SEKI1,2 , K. OGAWA1,2 , H. YAMAGUCHI1,2 , S. KAMIO1 , Y. FUJIWARA1 ,Y. KAWAMOTO1 , M. YOSHINUMA1 , T. KOBAYASHI1,2 , Y. TAKEMURA1 ,M. OSAKABE1,2 , M. ISOBE1,2 , M. YOKOYAMA1,2
12National Institute for Fusion Science, Toki, JapanThe Graduate University for Advanced Studies, SOKENDAI, Toki, Japan

Introduction The physics of the fast-ion transport is the one of the important issue for theachievement of the magnetically confined fusion reactor. This is because the plasma temperature is sustained by the kinetic energy transfer from fast-ions generated by external heatings andfusion born alpha particles. For this reason, the fast-ion transport mechanism has been investigatedin several magnetic confinement devices and simulations.
It is known that fast-ion confinement degrades with fast-ion-driven instabilities. The typicalfeature of this kind of loss is that the transport coeﬃcient depends on the fast-ion pressure, unlike other fast-ion loss mechanisms. This feature is explained that the amplitude of the instability
depends on the fast-ion pressure.
However, we observed the fast-ion confinement degradation depending on NB power even in
MHD quiescent plasmas in the Large Helical Device (LHD). In this paper, we estimate the confinement degradation quantitatively through the neutron measurement and the integrated simulation.
Experimental apparatus In this paper, two external heating systems, ECH and NBIs are used.
Plasmas are sustained by 3.8 MW ECH during each discharge. The magnetic configuration usedin this experiment are (Bt , Rax ) = (2.75 T,3.6m ), where Bt and Rax are the toroidal magnetic fieldstrength on the magnetic axis and the position of the magnetic axis with respect to the major radius,
respectively. The direction of the magnetic field is couter clockwise.
The LHD is equipped with five neutral beams as shown in Figure 1. In this work, we focus
on the tangential NB fast-ion. A positive ion-based perpendicular beam (NB#4) is used as a beamprobe for the charge exchange spectroscopy (CXS) measurements.
Experimentl Scenario Figure 2 shows the typical waveform used in these experiments. The
plasma discharges start at t = 3.3 s with 3.8 MW ECH. The first beam, NB#1, is injected from t =3.3 s to t = 5.3 s. The second beam, NB#3 (or NB#2), is injected at t = 4.3 s. In this paper, theperiod from t = 3.3 s to t = 4.3 s is called the “single beam phase” and the period from t = 4.3 s to t
= 5.3 s is called the “double beam phase”. In addition, the case with NB#1 and NB#3 is called the“co-injection case” and the case with NB#1 and NB#2 is called the “balance injection case”. NB#4
is injected repeatedly with a 20 ms pulse in each 200 ms from t = 3.32 s.
Experimental Results Figures 3 and 4 show the degradation of the fast-ion confinement due tothe increase in NB power. To consider the diﬀerence in the beam injection energy, the y-axis value
1

(NB#5) not used
SN183594
P nb [MW]6

port(NB#2) E~170keV
Pport ~2MW

δB/B0 [10 ] Β0=2.75Τ

Te(rho=0.5) [keV]
ne(rho=0.5) [e19 m-3 ](NB#1)
E~165keV
Pport ~2MW

neutron emission rate [e15/s](NB#3) E~170keV
Pport ~2MW

−3

(NB#4)
E~50keV
Pport ~6MW

Figure 1: Top view of NBI system in the LHD.
Typical beam injection energy and port-throughpower values are shown.

(a)

NB#1
NB#3
NB#4

4
2
0
3
(b)
2

Te
ne

1
0
1
(c)

0.8

183594 co-

0-20kHz

0.6
0.4
20-200kHz
0.2
0
1.6
(d)
1.2
0.8
0.4

double beam phase
NB#1 + (NB#3 or NB#2)

single beam phase
only NB#1

0

3.5

4

time [sec]4.5

5

Figure 2: Typical waveform used in these experiments (shot number (SN) 183594). (a) the NB
port-through power, (b) the electron temperatureand density, (c) the magnetic fluctuation, and (d)the neutron emission rate.
is Sn divided by the sum of the tangential NB port-through power multiplied by the beam-thermal
fusion cross-section and is normalized to that at t = 4.2 s.
Figure 3 shows that the neutron emission rate per beam power decreases by up to 20% after NB
superpose. On the other hand, the amplitude of the magnetic fluctuation does not change after NB
superpose as shown in Figure 2-(c). In addition, there are no significant variation in the electron
temperature and density and NB shine-through rate after NB superpose. For these reasons, thefast-ion confinement degradation causes the decrease of the neutron emission rate per NB power.
A similar reduction of the neutron emission rate per NB power appears in the balance injection casein the high-density discharges as shown in Figure 4. The neutron emission rate per NB power in the
double beam phase exceeds that in the single beam phase for the low-density discharges becausebalance injection cases (NB#1+NB#2)
Neutron emission rate per beam power
normalized to that at t=4.2s

Neutron emission rate per beam power
normalized to that at t=4.2s

co-injection cases (NB#1+NB#3)
1.4
1.2
1
0.8
0.6
0.4

nebar= 4e19
nebar= 3e19
nebar=2.5e19
nebar= 2e19

0.2
0
3.3

3.5

4

SN183595
SN183594
SN183591
SN183590

4.2
time [sec]nebar= 1.5e19 SN183586
nebar= 1e19 SN183585
nebar=0.75e19 SN183582

4.5

5

5.2

1.4
1.2
1
0.8
0.6
0.4

nebar= 3e19 SN183593
nebar= 2.5e19 SN183592
nebar= 2e19 SN183588

0.2
0
3.3

3.5

4

4.2
time [sec]nebar= 1.5e19 SN183587
nebar= 1e19 SN183584
nebar=0.75e19 SN183583

4.5

5

5.2

Figure 3: Time evolution of the neutron emis- Figure 4: Time evolution of the neutron emission rate per beam power for the co- injection sion rate per beam power for the balance injeccases.
tion cases.

2

the beam-beam fusion reaction is not negligible. From these result, we can guess that the fast-ion
confinement degradation also occurs in the balance injection cases and that the degradation may
not be visible due to the beam-beam fusion reaction in the low-density discharges.
Analysis model To quantitatively estimate the degradation depending on NB power, we formulate the following model. The time derivative of the neutron emission rate, S n , can be describedas
SnSn
PNB S n∂S n
= S − cl − amb −,(1)∂tτn τcxτano
where S is the sourve term of the neutron emission rate due to the beam injection, PNB is thedimensionless factor indicating the sum of the tangential NB port-through power normalized to 1
MW, τcln is the neutron decay time due to the classical slowing down, τamb
cx is the time constant of thecharge exchange loss, and τno is the time constant of the anomalous confinement degradation per 1
MW. The fourth term is the anomalous fast-ion loss, which is the focus of this paper. The source
term S can be roughly expressed as follows,
S ∼ Cnbulk
D PNB ,

(2)

where nbulk
D is the bulk-deuteron density and C is a constant factor related to the fusion cross-section
and the beam ion birth rate. The factor C can be viewed as a constant in a single discharge because
the beam injection energy, shine-through rate, plasma density, and temperature are constant in bothphases. By using values both in the single and double beam phases, we can obtain the following
expression about the time constant of the anomalous degradation.
τano =βαcl cl
PNB,d − nD,s
PNB,s )τamb( nD,dcx τn,s τn,dβ cl clβ amb cl
α cl clα amb clτ τ + nD,s
τcx τn,d − nD,dτn,s τn,d − nD,d
τcx τn,s
nD,s n,s n,d,(3)

where where subscripts “s” and “d” indicate the values in the single- and double-beam phases,α ≡ PNB,d /PNB,s , and β ≡ S n,d /S n,s , respectively. The eﬀective confinement time is given by()−1)( )−1 (
τano
amb −1eﬀ+.
(4)
τc,s = τcx
PNB,s
Simulation result with obtained τeﬀ
Figures 5 shows the time evolutions of the measured andcsimulated neutron emission rate for the co-injection case. Simulations with obtaind τeﬀ
c have a good
19
−3
19
−3
agreement in the range of 0.75 × 10 m < n¯e < 3 × 10 m . For this result, the obtained τeﬀ
c isvalid.
In contrast, eq. (3) can not be applied for the balance injection case because the beam-beamfusion is not ignorable. To apply the anomalous degradation for the balance injection case, thedependence of the anomalous degrdation on the electron density is estimated by the least-meansquare fitting approach. Figure 6 shows the dependece of Dano against the electron density, where
Dano = a2 /(5.8τano ) and a = 0.6 m. From the fitting, we obtain Dano = 0.049 × 10−19 ne .
By using the τeﬀ
c from the regression, the neutron emission rate simulation is performed forthe balance injection case. Figure 7 show the time evolution of the neutron emission rate for thebalance injection case. Simulation results have a good agreement even for the balance injection
case. From this result, the confinement degradation also occurs in the balance injection case.
3

1.4
1.2

SN183594_
ne =3.0*1019[m-3]τceff = 0.080s

1
0.8
0.6
0.4

measurementτceff = 0.15s, 0.080s

τceff = 0.15s

0.2
0
3.3

3.5

Dano :Diffusion coefficient [m2/s]

neutron emission rate [e15/s]

1.6

4

4.3
time [sec]

4.5

5

5.3

neutron emission rate [e15/s]

1.2
1
0.8

SN183582_
ne =0.75*1019[m-3 ]

0.6

0.4
Results in the co-injection casey=0.049x

0.3
0.2
0.1
0

0

0.5

0.4
measurementτceff = 0.32s, 0.21s

0.2
0
3.3

1

1.5
2
2.5
3
electron desity at ρpeak [1019 m-3]

3.5

4

Figure 6: The particle diﬀusion coeﬃcient due
to the anomalous loss against the electron denFigure 5: Time evolution of the measured and sity.
simulated neutron emission rate for the coinjection case.
3.5

4

4.3
time [sec]

4.5

5

5.3

1.6
neutron emission rate [e15/s]

neutron emission rate [e15/s]

1.6
1.4 SN183593
19 -3
1.2 ne (ρpeak )=2.6*10 [m ]
1
0.8
0.6
0.4

measurementτceff = 0.16s, 0.089s

0.2
0

3.5

4

4.5

5

1.4 SN183583
19 -3
1.2 n (ρpeak )=0.73*10 [m ]
1
0.8
0.6
0.4

measurementτceff = 0.32s, 0.22s

0.2
0

3.5

4

4.5

5

time [sec]time [sec]Figure 7: Time evolution of the measured and simulated neutron emission rate for the balance
injection case.
Conclusion We investigated tangential NB fast-ion confinement in MHD quiescent LHD plasmas. To investigate the degradation of the fast-ion confinement, we have performed a series ofexperiments using deuterium gas and deuterium neutral beams in the LHD. The experimental results show that the neutron emission rate per NB power decreases as the NB power is increased
in a single discharge. Because there is no significant change in the electron temperature, density,
and NB shine-through rate, the reduction in the neutron emission rate per NB power indicates thedegradation of the fast-ion confinement.
In this paper, we modeled the degradation to estimate the eﬀective confinement time based on
the measured data. The neutron emission rate simulation considering the confinement degradation
has a good agreement with the measured neutron emission rate. Therefore, the model, including
the confinement degradation depending on the NB power is valid.

4

Development of an adaptive predictive control system using
data assimilation at LHD
Y. Morishita1, S. Murakami1, N. Kenmochi2,4, H. Funaba2, M. Yokoyama2,4, G. Ueno3,4,5, M.
Osakabe2,4
1Department of Nuclear Engineering, Kyoto University, Nishikyo, Kyoto, 615-8540, Japan
2 National Institute for Fusion Science, National Institutes of Natural Science, Toki 509-5292, Japan
3The Institute of Statistical Mathematics, Research Organization of Information and Systems, Tokyo 190-8562,

Japan
4The Graduate University for Advanced Studies, SOKENDAI, Toki 509-5292 and Tokyo 190-8562, Japan
5The Joint Support-Center for Data Science Research, Tokyo 190-0014, Japan

1. Introduction
The operation of future fusion reactors requires a system that predicts and controls
fusion plasma behavior under large uncertainties. To develop a numerical system that cancontrol the behavior of fusion plasmas with high accuracy, we are developing a data
assimilation system, ASTI [1], based on the integrated transport simulation code, TASK3D[2]. Data assimilation is a statistical method to optimize numerical simulation using
observations. To realize a control system that can adapt the predictive model to the actual
plasma and estimate the optimal control input based on the prediction, we have extended thedata assimilation framework to include control procedures (DACS [3]) and implemented anadaptive predictive control system at the Large Helical Device (LHD) [4]. In this paper, weshow the result of LHD experiments to control the central electron temperature and the issues
to more advanced control experiments.

2. Data assimilation-based control system in LHD
ASTI is a data assimilation system which employs the integrated transport simulation
code, TASK3D, as the system model. We have extended ASTI to a control system based on
the DACS framework. ASTI assimilates the experimental time series data into the systemmodel (TASK3D) to get the model behavior close to the actual behavior and estimate thecontrol input for the predetermined target state using the optimized system model. ASTI canprovide an adaptive model predictive control for fusion plasma.
As a first step to demonstrate the control performance of the system, we performed
LHD experiments to control the central electron temperature. ASTI adjust the ECH total
power to bring the central electron temperature to the target temperature. Four gyrotrons areused for the ECH; thus, ASTI controls the ECH power by selecting a subset of these four
gyrotrons. ASTI sends the estimated control signals to the gyrotrons every 0.1s andassimilates the observed radial profiles of electron temperature and density every 0.3s. The
radial profiles of the electron density and temperature obtained by the real-time Thomson
scattering measurement system are assimilated into the system model to reduce theuncertainty in the electron and ion thermal diffusivities and improve the prediction
performance of the system model.

3. Results and discussion
Figure 1 (a) shows a control result of electron temperature at the plasma center for thetarget temperature of 4 keV, and Fig. 1 (b) shows the adjusted ECH power. We can see thatthe electron temperature approaches the target temperature at 3.9 s. The control error(differences between the observed temperature and the target temperature) increases to 1 keVin the transient section (2.1-3.9 s) due to the imperfection of the system model. However, theerror decreases to about 0.3 keV in the steady state section (3.9 s-). This is because thethermal diffusivity models have been adjusted to bridge the gap between the TASK3D
simulation and the actual LHD plasma in the adaptation step.

Fig. 1 Results of a control experiment (#186502). (a) Control result of the electron
temperature at the plasma center. (b) Adjusted ECH power by ASTI. The blue hatchingareas represent the one standard deviation of the predicted distribution.
We must address the following issues for more advanced experiments in the future.
• To prevent a reduction in the control performance due to delays in adaptation, we require acontrol algorithm that reflects observed information into the prediction and control
estimation more quickly.
• The probability distribution of model parameters should be adjusted by considering thephysical characteristics and control accuracy to avoid control instability.
• Reduction of the control error in the transient phase makes it necessary to improve theassumed transport model.

Referencesfi[1] Y. Morishita, et al., Computer Physics Communications 274, 108287 (2022).
[2] S. Murakami, et al., Plasma Phys. Control. Fusion 57, 119601 (2015).
[4] Y. Morishita, et al., J. Computational Science 72, 102079 (2023).
[5] Y. Morishita, et al., Scienti c Reports 14, 137 (2024).

Study of model-based control fortokamak plasma profile
M. Sato1, M. Kishida1,2, S. Sakakibara1,3, *M. Yokoyama1,3
1
SOKENDAI, 2National Institute of Informatics, 3National Institute for Fusion Science

1. Introduction
Controlling the plasma profile remains a significant challenge in achieving afusion reactor. To control highly autonomous complex dynamics, such as burning
plasmas, conventional PID control mechanisms are inadequate. Therefore, we are
seeking ways to control the plasma profile from the viewpoint of control theory andaim to develop more advanced control methods.

2. What is control theory?
Control theory is the study of systems that have inputs and outputs, and it involves
modeling, analyzing, and finally designing controllers.
The first step in controller design is modeling. There are three main methods ofmodeling: first-principles modeling, black box modeling, which only uses therelationship between inputs and outputs, and grey box modeling. System models
are typically defined using state equations. It is expressed as follows:𝒙̇ = 𝒇(𝒙, 𝒖) (1)𝒚 = 𝒈(𝒙, 𝒖) (2)where 𝑥 ∈ ℝ! is state vector, 𝑢 ∈ ℝ" is input vector and 𝑦 ∈ ℝ# is output
vector. There are different ways to express the state equation, including linear
time-invariant systems (LTIs), nonlinear systems, etc.
The second step is system analysis. System stability analysis is an essential aspectthat needs to be given special attention. This is because the major purpose ofcontrol is to stabilize the system.
The last step is “tools” selection. Control theory encompasses various tools, suchas Robust control theory which considers model uncertainty, and Optimal control
theory which deals with constraints. When designing controllers, several tools canbe used for various purposes.

3. Plasma profile control model
Plasma profile control models have been developed globally using different
approaches. For example, a 1D diffusion-like parabolic equation for the radial
profile of the plasma poloidal flux or current density. Readers could refer to Witrant
et al. [2007]. This is a control model known as a distributed parameter system incontrol theory. It requires advanced mathematical understanding. Another example
is a port-Hamiltonian system for plasma profile which is referred to in Vu [2012].
A port-Hamiltonian system is a framework used in physics and control theory tomodel and analyze complex dynamical systems. Various profile control
experiments are being conducted around the world using these control models.

4. Prospects for future studies
To control complex systems, advanced control logic is required; particularly, theport-Hamiltonian formulation of burning plasma will be studied intensively for
future studies. In recent years, there have been many data science approaches,
including AI. However, it is always necessary to have a control theory framework
when discussing system stability and other issues. In fact, since various instabilities
exist in plasmas, such as MHD instability, it is crucial to control these instabilities.
A comprehensive understanding of various factors and plasma control is essential
to realize future fusion reactors. Model-based control is a promising solution.

5. References[1] Witrant, E., Joffrin, E., Brémond, S., Giruzzi, G., Mazon, D.,Barana, O., andMoreau, P. (2007). A control-oriented model of the current profile in tokamak plasma.
Plasma Phys. Control.Fusion, 49, 1075–1105.
[2] Vu, NM., Lefevre, L., Maschke, B. (2012). Port-Hamiltonian formulation forsystems of conservation laws: application to plasma dynamics in Tokamak reactors.
IFAC Proceedings Volumes Volume 45, Issue 19, 2012, Pages 108-113

Prediction of radiative collapse by analyzing videos ofplasma discharge with CNN in LHD
Yuya Suzuki1, Mamoru Shoji1,2, Naoki Kenmochi1,2, and Masayuki Yokoyama1,2
1

SOKENDAI, 2National Institutes of Natural Sciences, National Institute for Fusion
Science

1. Introduction
The unexpected termination of plasma discharges has been a critical issue fornuclear fusion research. One of the termination events observed in the Large Helical
Device (LHD) is named “radiative collapse”. This event occurs in the plasma edge
region and is caused by impurity influx, which cools down the plasma temperature
with strong luminescent light and prevents the sustainment of high-density plasmas[1]. To avoid this event, a new method that predicts the radiative collapse from theimages of plasma discharges is proposed using Convolutional Neural Networks(CNN) which is one of the machine learning models.
2. Dataset for training model
To train a model on features of radiative collapse, a dataset was created fromdischarges with radiative collapse. The search was conducted using experimental datawith a magnetic axis of 3.6 m, a magnetic field of 2.75 T, and shot numbers
ranging from 157000 to 162000 with hydrogen as the fuel. As a result, 41 discharges
were selected. From plasma discharges video of these discharges, images of stable
and collapse state were extracted and labeled. These images were used for dataset.
3. Accuracy evaluation of trained determination models
A CNN model called EfficientNet[2] was used to determine the discharge statefrom the images. The training results of the model for each camera are shown infigure 1. Except for a few cameras, the accuracy during training, accuracy on theevaluation dataset and F-score were all above 86%. This result shows that themodels in those cameras were successfully trained and a determination model withhigh accuracy was obtained.

Fig 1. Results of trained models at
various evaluation parameters. Outerlabels show camera’s name.
Learning : Accuracy on the training dataset
Demo : Accuracy on the evaluation dataset
F score : Accuracy considering data biasMiss : Misjudged images in collapse state
False : Misjudged images in stable state.

4. Performance of determination models as predictor
The determination model was applied to a plasma discharge video as a predictor.
The results are shown in Figure 2. The blue hatched area is the area defined asstable, so models with a collapse judgement from this section are excluded. As aresult, models such as 10-O and 2-O were able to predict about 2-4 frames (66-132
milliseconds) before the radiative collapse.

Fig 2. Average outputof the discharge state
determination modeland time to onset ofradiative collapse at
the point of output.

5. Analysis of the focus area of the determination models
The determination basis of the model was visualized by a method called
GradCAM[3]. The results are shown in Figure 3. The results show that there were
two main determination patterns of the model: one was the emission of plasma near
the divertor tiles, and the other was the transparency of the plasma. As these features
are areas that have been mentioned in previous studies [1], it was found that themodel focused on significant areas in the image to determine.

Fig 3. Heat map showing the area of interest for the discharge state determination model by
GradCAM.

6. Conclusion
A discharge state discrimination model could be created by using camera images
accumulated from previous plasma discharge experiments at LHD for machinelearning. For the cameras with good accuracy, the determination accuracy of theimages was about 86-92% and the prediction performance was about 2-4 frames,
which was sufficient for predictive control. These models were also found to
determine from areas in the image that captured phenomena such as thosementioned in previous studies. The results demonstrate the effectiveness ofmachine-learning-based controllers and data-driven research in plasma and fusion
research.
7. Reference[1] B.J. Peterson, et al., Plasma Fusion Res. 1, 045, (2006).
[2] Mingxing Tan and Quoc V.Le, arXiv:1905.11946 (2019)[3] Selvaraju R.R., et al., arXiv:1610.02391 (2016)

Introduction of hydrogen isotope eﬀects to a turbulent
transport model DeKANIS and its application to integrated
simulations
E. Narita1 , M. Honda1 , M. Nakata2,3 , N. Aiba4 , N. Hayashi4 , T. Nakayama3 , andM. Yoshida4
1

Graduate School of Engineering, Kyoto University, Kyoto, Japan
2
National Institute for Fusion Science, Gifu, Japan
3
The Graduate University for Advanced Studies, Gifu, Japan
4
National Institutes for Quantum Science and Technology, Ibaraki, Japan

1

Introduction

Hydrogen isotope eﬀects are crucial in predicting burning plasma performance with the knowledge of physics in deuterium and hydrogen plasmas. Experimental observations show positive
isotope eﬀects in several devices (see e.g. [1,2]), and several factors have been reported withgyrokinetic codes to explain the experiments (see e.g. [3,4]). In this study, we have introducedhydrogen isotope eﬀects into a turbulent transport model DeKANIS [5] and performed temperature and density proﬁle predictions with the modiﬁed model.

2

How to introduce the hydrogen isotope eﬀects into DeKANISDeKANIScalculates the )electron particle( and heat ﬂuxes with) the following expressions:
(
R
R
R
R
Γ̄e = D̄ Ln + CT LT + CP and Q̄e = χ̄e CN Ln + LT + CHP . Here, non-dimensional D̄,eeeeχ̄, R/Ln and R/LT are the particle and heat diﬀusivities and the normalized density and temperature gradients, respectively. The ﬂuxes consist of the diagonal and oﬀ-diagonal terms,
which are quantitatively determined by estimating the oﬀ-diagonal term coeﬃcients (CT,P,N,HP )and the diﬀusivities (D̄ and χ̄e ). DeKANIS also predicts the ion heat ﬂux, which is given asχ̄i,eﬀχ̄e,eﬀ LRT nnei TTei . Here, χ̄eﬀ is the normalized eﬀective heat diﬀusivity. Once χ̄e,eﬀ isQ̄i = χ̄e,eﬀ
i

obtained with Q̄e , Q̄i can be calculated by determining χ̄i,eﬀ /χ̄e,eﬀ . To forecast Γ̄e , Q̄e and Q̄i ,DeKANIS utilizes a neural network (NN) model and a saturation model. The NN model istrained to evaluate CT,P,N,HP and χ̄i,eﬀ /χ̄e,eﬀ , whose training datasets are prepared by the linear
gyrokinetic calculations with GKW [6] at the normalized wavenumber k̄θ where the normalized
linear growth rate takes the maximum value, γ̄max , of the ion temperature gradient (ITG) modeor the ITG/trapped electron mode (TEM). The NN model also outputs the corresponding k̄θ,γ̄max
and γ̄max , which are passed to the saturation model to determine D̄. The remaining coeﬃcient,χ̄e , is calculated with the relationship between D̄ and χ̄e : χ̄e = D̄(CT + 1.5)/CN .
The previous training datasets of the NN models were given only by the calculations assuming
deuterium plasmas. To consider the changes in the output parameters due to the main ion mass
mi , we have introduced mi as an input parameter and added the calculations assuming hydrogenand tritium plasmas to the training datasets. The saturation model has also been modiﬁed

Figure 1. Wavenumber spectra of quasilinear diﬀusivities for hydrogen, deuterium and tritium
plasmas. The circles denote the quasilinear diﬀusivities at the wavenumber where the linear
growth rate takes the maximum value. (b) The ratio of the wavenumber integral of the quasilinear diﬀusivity to the quasilinear diﬀusivity at the wavenumber where the linear growth ratetakes the maximum value as a function of the wavenumber. (c) Comparison of the wavenumber
integral of the quasilinear diﬀusivity and its model. (d) Comparison of the semi-empiricallyestimated particle diﬀusivity and its model for hydrogen and deuterium plasmas.
founded on the hydrogen isotope eﬀects given by linear gyrokinetic calculations. The details of
the new saturation model is shown in the next section.

3

Modiﬁcation of the saturation model

The new saturation model includes the two kinds of hydrogen isotope eﬀects. The ﬁrst one
is the diﬀerence in the suppression of the TEM due to collisions depending on mi [4]. When
we perform linear gyrokinetic calculations and compute the quasilinear diﬀusivity γ̄/k̄θ2 as afunction of k̄θ like in Fig. 1(a), the wavenumber integral of γ̄/k̄θ2 can be estimated as a proxy of
a diﬀusivity given
calculation. Here, the wavenumber integral of γ̄/k̄θ2 is deﬁned[∑ (by a nonlinear])2
as ⟨γ̄/k̄θ2 ⟩ =k̄θ γ̄/k̄θ ∆k̄θ /k̄θ,lim , where ∆k̄θ is the interval of k̄θ used for the calculationand k̄θ,lim is the upper limit of the integration. In this study, k̄θ and γ̄, are normalized withthe deuterium unit. Since the TEM is more suppressed for the heavier ion mass case due
to collisions [4], ⟨γ̄/k̄θ2 ⟩ decreases with the increase in mi as found in Fig. 1(a). The second
hydrogen isotope eﬀect is the change in the residual zonal ﬂow level. When we perform linear
gyrokinetic calculations, a reduction in the residual zonal ﬂow level with an increase in thecollision frequency is found. This tendency can be expected by [7], where the suppression of
2 τ /ν . Here, ω =turbulence by zonal ﬂows considering collisional decay is given as νf2 /ωf4 + KRH
c d
f
( 2)
vth,i /R, νf = ωf exp −q , KRH is the linear zonal ﬂow response function without collisions [8],τc is the correlation time and νd = 0.67/ (τii ϵ) with the ion thermal velocity vth,i , the major
radius R, the safety factor q, the inverse aspect ratio ϵ and ion-ion collision time τii . A square2 τ /ν , which is deﬁned as characteristic time of the residual zonal ﬂow, τ ,root of νf2 /ωf4 + KRH
c dr
in this study, is found to be proportional to m0.5.
i
The second hydrogen isotope eﬀect can be included by introducing τr to the saturation model,
but the inclusion of the ﬁrst one requires a model that expresses ⟨γ̄/k̄θ2 ⟩ with k̄θ,γ̄max and γ̄max .
2k̄θ,γ̄max tends to be lower for the heavier ion mass case, and γ̄max /k̄θ,γ̄
often does not reﬂect themax
2tendency of ⟨γ̄/k̄θ ⟩ as shown with the circles in Fig. 1(a). The relationship between ⟨γ̄/k̄θ2 ⟩ and2γ̄max /k̄θ,γ̄
was investigated by performing linear calculations with diﬀerent mi for a variety ofmax
2
plasma parameters with k̄θ,lim = 1.2, and it was found that ⟨γ̄/k̄θ2 ⟩/(γ̄max /k̄θ,γ̄) increases withmaxk̄θ,γ̄max as shown in Fig.
1(b). This
with the following functional formula:() (trend can be expressed)α22⟨γ̄/k̄θ ⟩model ∼ coef γ̄max /k̄θ,γ̄max k̄θ,γ̄max /k̄θ,lim . coef and α were determined by regression
analysis utilizing the calculation results shown in Fig. 1(b) as coef = 3.01 and α = 1.33, andthe obtained model can express ⟨γ̄/k̄θ2 ⟩ for all hydrogen isotope species (Fig. 1(c)).

Figure 2. Radial proﬁles of the electron and ion temperatures and the electron density predictedby the transport simulation with TOPICS for ITER PFPO-1 relevant plasmas, with the change
in the main ion mass.
[

]1.20

0.752 0.1τ̄ ⟨γ̄/k̄ 2 ⟩0.5
The new saturation model was obtained as D̄ = 1.23⟨γ̄/k̄θ2 ⟩model. Therθ modelconstant 0.1 in the zonal ﬂow part adjusts its order. The constant 1.23 and the exponents were
determined by regression analysis to reproduce the semi-empirically estimated particle diﬀusivity D̄exp . Although D̄exp was estimated based on electron heat ﬂuxes of JT-60U deuterium
plasmas, the obtained saturation model can express D̄exp for both the deuterium cases used forthe regression analysis and the hydrogen cases introduced as test cases as shown in Fig. 1(d).
Therefore, the saturation model could be used also for tritium plasmas.

4Integrated simulations including the hydrogen isotope eﬀects

The modiﬁed DeKANIS was implemented in the integrated code TOPICS [9] and temperatureand density predictions were performed for plasmas envisioned in the ITER pre-fusion power
operation 1 (PFPO-1) phase with three cases having diﬀerent main ion species: hydrogen,deuterium and tritium, as shown in Fig. 2. Here, the magnetic equilibrium and the pedestal
height are ﬁxed, the plasma current is 5.0 MA, the toroidal magnetic ﬁeld is 1.8 T and theelectron cyclotron heating power is 30 MW. The electron temperature and density are higher
for the heavier ion mass case, but the ion temperature is higher for the lighter ion mass case.
The change in the electron proﬁles should be caused by the reduction in the diﬀusivities with
the increase in mi , and that in the ion one may be due to the collisional equipartition, which islarger for the lighter ion mass case. The thermal energy conﬁnement time is 2.15, 2.34 and 2.42
for the hydrogen, deuterium and tritium cases, respectively. The predicted dependence of thethermal energy conﬁnement time on mi is weaker than that of the scaling the IPB98(y,2) scaling.
It may be due to the same pedestal height for the three cases. We will develop a pedestal model
including hydrogen isotope eﬀects in the future.

References[1] H. Urano et al., Nucl. Fusion 53 083003 (2013).
[2] C.F. Maggi et al., Plasma Phys. Control. Fusion 60 014045 (2018).
[3] J. Garcia et al., Plasma Phys. Control. Fusion 64 054001 (2022).
[4] M. Nakata et al., Phys. Rev. Lett. 118 165002 (2017).
[5] E. Narita et al., Contrib. Plasma Phys. e202200152 (2023).
[6] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[7] F.L. Hinton and M.H. Rosenbluth, Plasma Phys. Control. Fusion 41 A653 (1999).
[8] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80 724 (1998).
[9] N. Hayashi and JT-60 Team, Phys. Plasmas 17 056112 (2010).

Research of turbulent transport due to dissipative trapped electronmode in tokamak plasmas
S. Toda1,2 , M. Nunami1,3 , and N. Kasuya4
1 National Institute for Fusion Science
2 Graduate Institute for Advanced Studies, SOKENDAI
3 Graduate School of Science, Nagoya University
4 Research Institute for Applied Mechanics, Kyushu University

Turbulent transport is one of the most important subjects in research of toroidal plasmas. The
purpose of this study is to quantify turbulent transport. Microscopic plasma instabilities arestudied using electromagnetic gyrokinetic simulations. The trapped electron mode (TEM) where
normalized electron collision frequency is less than unity 1 is called dissipative-TEM in a high collisional regime. Here, this mode is driven by dissipation from the electron collision and ions do
not need to be in the banana regime. This mode is known to be driven unstable in the presence ofdensity and/or temperature gradients. Plasma experiments at the PLATO device at the ResearchInstitute for Applied Mechanics, Kyushu University have been started. Transport simulations using integrated code are carried out on the PLATO tokamak to forecast plasma performance andplasma proﬁles in the PLATO are predicted using an empirical model in the integrated code, the
TASK code. In these plasmas, the electron collision frequency is typically less than the electron
bounce collision frequency and the ion collision frequency is larger than the ion bounce collision
frequency. Therefore, these plasmas are in a high collisional regime and the d-TEM is predicted
to be unstable. We focus on the d-TEM and the ITG mode by using local gyrokinetic simulations.
The GKV code2 is used for local ﬂux tube gyrokinetic simulation, using the Sugama collision model
operator3 or the Lenard-Bernstein collision model operator. The Sugama collision operator wasused to study the linear response of the zonal ﬂows with a collisional eﬀect 4 . In these analysis
results, the residue level of the zonal ﬂows using the Sugama collision operator was found to beclose to the analytically predicted one. On the other hand, the linear response of the zonal ﬂows
using the Lenard-Bernstein collision operator shows over-damping. Therefore, the Sugama collision model operator is used for the local gyrokinetic simulation in this study.
The plasma proﬁles predicted by the TASK code in the PLATO are used. The linear simulation
results using the Sugama (S) and Lenard-Bernstein (LB) collision operators are compared. The
linear growth rate is studied at kx ρe = 0.0 and in 0.0 ≤ ky ρe ≤ 0.06 region, where kx and ky are
the radial and poloidal wavenumbers, respectively. The dependences of linear growth rate γ, andreal angular frequencies ω, on the normalized poloidal wavenumber ky ρe , are studied. The linear
growth rate γ, and the real angular frequencies ω, as the function of ky ρe are also investigated.
Using the S collision operator, the instabilities predicted in this study are driven by the d-TEM
in the regions 0.005 ≤ ky ρe ≤ 0.020 at ρ = 0.47 and 0.005 ≤ ky ρe < 0.040 at ρ = 0.65. The real
angular frequency changes from the electron to ion diamagnetic drift motion direction at ρ = 0.47
and ρ = 0.65 when the poloidal wavenumber increases. The so-called d-TEM is unstable, where
the electron collision frequency is less than the electron bounce frequency and the ion collision
frequency is larger than the ion bounce frequency1 . This d-TEM is driven by the density gradientin this study. The d-TEM is predicted to be unstable at ρ = 0.81 in the all-poloidal wavelength
region studied here. The linear growth rate for the d-TEM at ρ = 0.81 becomes larger than that at
ρ = 0.47 because of the increase of R/Ln , therefore, the d-TEM becomes more unstable in the outerradial region. On the other hand, the ITG mode is excited in the regions 0.020 < ky ρe < 0.060 at
ρ = 0.47 and 0.040 ≤ ky ρe < 0.060 at ρ = 0.65. The ITG is stable at ρ = 0.81 because R/LTi issmaller than that at ρ = 0.47 or ρ = 0.65. The rapid changes from the d-TEM to the ITG mode are
obtained in the dependence of the real angular frequency on the poloidal wavenumber, specially atρ = 0.65. Using the LB collision operator, only the d-TEM is unstable in the 0.005 < ky ρe < 0.030

region at all radial points. The diﬀerences between the results using S and LB collision operators
are found. These are because the ﬁeld particle part is included with momentum conservation andthe collisionality term depends on the velocity in the S collision model operator3,4 .
Next, we study the linear response of the zonal ﬂows using the S and LB model collision operators. The zonal ﬂow decay time5 is investigated by the linear gyrokinetic simulation results. We
examine the linear response of the zonal ﬂows using S and LB collision
operators.
⟨⟩ ⟨ Here, the linear⟩
response function of zonal ﬂows Rk̃x (t) is deﬁned as Rk̃x (t) = ϕ̃k̃x ,k̃y =0 (t) / ϕ̃k̃x ,k̃y =0 (t = 0) ,where ϕ̃(= ϕ/(Te ρe /(eR))) is the electrostatic potential ﬂuctuation. We set k̃x = 0.01. The residue
levels at ρ = 0.47, ρ = 0.65 and ρ = 0.81, which is deﬁned as the value of Rk̃x (t = 2500) here, arestudied. The residue level using the∫S collision operator is larger than that using the LB one. The
zonal ﬂow decay time is deﬁned as 02500 Rk̃x (t)dt in this article. The zonal ﬂow decay time is alsostudied. The zonal ﬂow decay time using the S collision operator is found to be twice as large asthat using the LB collision operator. The linear simulation results suggest that the eﬀect of thezonal ﬂows on the turbulent transport using the S collision operator in the nonlinear simulation
results is larger than that using the LB collision operator.
A nonlinear gyrokinetic simulation is carried out. A nonlinearly saturated state for plasma turbulence is obtained. Nonlinear gyrokinetic analysis
results⟨are studied(⟩ )for the time evolutions of the2squared turbulent potential ﬂuctuation, T = Σk̃x ,k̃y ̸=0 |ϕ̃k̃x ,k̃y | /2 and the squared zonal ﬂow
(

⟨

⟩

)

potential, Z = Σk̃x ,k̃y =0 |ϕ̃k̃x ,k̃y |2 /2 . Time evolutions for the squared potential ﬂuctuations using the S and LB collision model operators are compared at ρ = 0.47, ρ = 0.65 and ρ = 0.81. To
study the eﬀect of the zonal ﬂows, the values of T̄ and Z̄ at ρ = 0.47, 0.65 and 0.81 are examined.
The symbol¯represents the averaged values of T or Z of time t normalized by R/vte from 3000 to
5000. The eﬀect of the zonal ﬂows on the turbulence by the S collision operator is predicted to
be much larger than that by the LB one. The values of Z̄/T̄ at ρ = 0.47, 0.65 and 0.81 using the
S and LB collision model operators are studied. The values of Z̄/T̄ using the S collision model
operator are larger than those using the LB collision model operator. Next, the transport level isevaluated by the nonlinear simulation results. The averaged values over the time interval after the
nonlinear saturation for Qe , χe , Qi , χi , Γ and D using the S and LB collision model operators are
obtained, where Qj is the energy ﬂux and χj is the energy diﬀusivity for the species j. Here, ΓGBGBis particle ﬂux and D is particle diﬀusivity. The normalized values of Qe /QGBe , χe /χe , Qi /Qe ,GBGBGBχi /χe , Γ/Γ and D/χe using the S and LB collision model operators are also obtained. Energy
and particle ﬂuxes at ρ = 0.47, 0.65 and 0.81 are studied. The energy and particle diﬀusivities at
ρ = 0.47, 0.65 and 0.81 are indicated. The energy and particle ﬂuxes normalized by the gyro-Bohm
unit at ρ = 0.47, 0.65 and 0.81 are also studied. The energy and particle diﬀusivities normalized
by the gyro-Bohm unit at ρ = 0.47, 0.65 and 0.81 are also investigated. The transport levels using
the S collision model operator are found to be smaller than those using the LB one except in thecase of electron energy transport at ρ = 0.81 due to the eﬀect of the zonal ﬂows.

REFERENCES
1

B. B. Kadomtsev and O. P. Pogutse, Reviews of Plasma Physics, edited by M. Leontovich, Vol. 5(consultants Bureau, New York, 1970).
2
T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
3
H. Sugama, T. H. Watanabe, and M. Nunami, Phys. Plasmas 16, 112503 (2009).
4
M. Nunami, M. Nakata, T. H. Watanabe, and H. Sugama, Plasma and Fusion Research 10,1403058 (2015).
5
S. Ferrando-Margalet, H. Sugama, and T.-H. Watanabe, Phys. Plasmas 14, 122505 (2007).

統合乱流計測シミュレータの開発とそのトーラス装置計測への応用
Development of Integrated Turbulence Diagnostic Simulator and Its
Application to Torus Device Measurements
N. Kasuya1,2, R. Yoshihara2, T. Ido1,2, M. Nunami3, K. Tanaka1,3, S. Satake3, K. Fujita4,
M. Yagi5
1
RIAM, Kyushu University, Kasuga, Fukuoka, Japan
2
IGSES, Kyushu University, Kasuga, Fukuoka, Japan
3
National Institute for Fusion Science, Toki, Gifu, Japan
4
Gradutae School of Science, Nagoya University, Nagoya, Aichi, Japan
5
National Institute for Quantum and Radiological Science and Technology, Rokkasho, Aomori,JapanE-mail: kasuya@riam.kyushu-u.ac.jp
1. Introduction
Plasma turbulence has large effects on particle and heat transport in magnetically-confined
fusion plasmas, and it is important to clarify its 3-D feature in real magnetic geometries intorus experiments, because various global and meso-scale structures are formed. Synthetic
diagnostics use simulation data to represent experimental measurements by simulating
measurement signals [1-5]. Numerical simulations can show global features of mean profile
modifications and fluctuation excitations, and quantitative comparison between experiments
and simulations is possible. We have been carrying out some numerical diagnostics forturbulence and structural formation analyses [6,7]. On the other hand, several kinds ofsimulations use different models, depending on the simulation targets, so comparison between
models is also important to capture various characteristics in magnetized plasmas. In addition,
combination of multiple simulation results can give complex plasma characteristics including
several physical mechanisms. Therefore, we newly develop these analysis routines into an
integrated numerical diagnostic platform for multiple simulation data [8].
2. Integrated turbulence diagnostic simulator
The integrated turbulence diagnostic simulator (iTDS) is a framework for making global
analyses in real magnetic geometries of fusion devices [8,9]. iTDS connects global
simulations and experimental observations. Various analysis routines for several experimental
diagnostics are developed on this platform. The analysis consists of the following processes;i) 3-D simulation of plasma instabilities and transportii) Plasma device parameter set and adjustment of the data format (if necessary)iii) Simulation to obtain signals by experimental diagnostics
In process i), several codes can be used to obtain 3-D plasma profiles and time evolution offluctuations. Global codes are used for calculations in torus plasma geometries. In general,
particular models do not always include all the necessary quantities, so lacked information is
supplemented with other calculation or assumption, if necessary in process ii). In addition,
spatial coordinates in simulation codes may be different with ones used for numericalmeasurement, so data transfer is necessary in that case. Data mapping to the common format
is often used, but here we introduce the object-oriented method in accordance with each data
format to gain the maximum information. The data transfer is limited to the minimum. Thenthe numerical measurements are carried out in process iii). The basic component to capture

3-D dynamics from simulation data is extracting the local value at the arbitrary position fromdiscreet data by interpolation. The combination of this basic routine can give multi-point
observation along a line of the sight of each diagnostic. Here the focus is mainly set on ageometric effect from magnetic configurations. Applications of iTDS is explained in thefollowings to show mean profile and fluctuation measurements in torus experimental devices.
3. Steady profile measurement
HIBP is powerful diagnostic to give information of several kinds of quantities [10]. Single
trajectory of the detected ion can be selected with given injection and detection conditions,which determines the observation position. Plasma fluctuations and electromagnetic fields
inside the plasma can be measured simultaneously without disturbing the plasma. Calculationof 3-D trajectory of heavy ions, including non-circular plasma shape and experimental
arrangement of diagnostics, is carried out by synthetic HIBP (sHIBP) in iTDS for LHD. The
following set of equations of motion in cylindrical coordinates (R, , z) is used for the
calculations;  R  2  q   q RB  zB Rzm R m1
q 1  q 1
 
 R  RB  2 R 

 zB
z
R
m R 2  m R
q  q 1 zRB  RBRm z m R

 









(1a )
(1b)
(1c)

Mean radial 0 and 3-D perturbation 1 profiles of the potential are included in the 3-D
magnetic configurations with (BR, B, Bz). The 3-D variation of the potential is calculated by
global drift kinetic code FORTEC-3D, which is given by the neoclassical process [11]. This
3-D variation is predicted to affect the impurity neoclassical transport. An example of thecalculated trajectories of the injected heavy ions are shown in Fig. 1 with their ionization
points from charge number +1 to +2, which are the observation points for HIBP. For the
calculations, it is necessary to use magnetic fields for the whole regions inside and outside ofthe last closed flux surfaces (LCFS). The magnetic fields generated by external coils (vacuum
fields) are modified to be equilibrium fields in the plasma discharges. FORTEC 3-D solves
the set of equations by using magnetic fields from equilibrium code VMEC [12] only inside

FIG. 1. HIBP trajectory from the injection point to the detection point in the vacuum vessel.
Three cases are indicated with usage of (i) vacuum magnetic fields for all regions, (ii)equilibrium magnetic fields inside of LCFSs, and (iii) equilibrium magnetic fields and 3-D
potential variations. The circles indicate the ionization points of injected ions, whose spatialdifferences are less than 1cm as in the magnification.

of the LCFSs, as is usual for plasma simulations. On the other hand, the modification of themagnetic fields is small outside of the LCFSs, so it is convenient to use the vacuum fields tocover wider region in the vacuum vessel. In this sHIBP, we combine to use the equilibrium
fields inside of the LCFSs and the vacuum fields outside of the LCFSs. This combination iseasy to be implemented by using the module structure of iTDS. Figure 1 shows thecomparison of the cases with the combined fields and only with the vacuum fields. The
difference is rather small, considering the resolution of experimental detectors. The
comparison between the cases with and without potential variation 1 is also included in Fig.
1 to confirm little effect of the perturbed electric fields on the heavy ion trajectories. This isbecause the injected ion energy (several MeV) is much larger than the plasma potential(several 10 eV for the maximum). In this way, the effect of the 3-D variations is evaluated tofind out the observation points of HIBP. Interactive evaluations of density and magnetic
fluctuations are also possible, which are important to identify errors of experimental
observations.
In sHIBP, variation of the beam energy (velocity of the heavy ions) and attenuation rate of
the beam intensity are evaluated to give information of the plasma potential and density along
the HIBP trajectory, as in Fig. 2. The magnetic fields (b), electrostatic potential (c) andelectric fields (d) (spatial derivatives of the potential) are input parameters for the sHIBP, inwhich the set of the equations of motion is solved to calculate the trajectories. The connectionof the vacuum and equilibrium fields at the LCFSs gives little discontinuity in the magnetic
fields as in Fig. 2(b). The moment of the ionization from charge number +1 to +2 is indicatedby the dashed line in Fig. 2. Figure 2(e) shows change of the heavy ion velocity along the
trajectory, corresponding to the plasma potential profile in Fig. 2(b). The sudden change of theeffective potential at the ionization point gives difference of the energies between the injected
and detected heavy ions as in Fig. 2(e), which gives observation of the potential value at theionization point. Noted that the non-smooth lines in Fig. 2(d) is due to usage of data on themeshes with finite intervals, which cause numerical errors for calculations. The 0 variation
is 10-3 orders of magnitude smaller than the injected beam energy, so is observable inexperiments, but the 1 variation from FORTEC-3D simulation is 10-6 orders of magnitude
smaller than the injected beam energy, so is too small for detection. The detected beam
intensity is given as
I S  I p Ap As ps vrvb

n  rs  l ,(2)

where I is the beam intensity, vb is the beam velocity, vr is the relative velocity, ps is thecollisional cross-section of the ionization from the primary beam to the secondary beam, n(rs)is the density at the ionization point, l is the effective reaction path length of the ionization
process, and the lower suffixes i = p and s denote the primary and secondary beams,
respectively. Beam attenuation rate Ai is given with
 i vr
Ai  exp   
ndl  ,
Li
vb



(3)

where Li is the path length along the trajectory. Figure 2(a) shows attenuation rates of the
primary and secondary beams (ln Ai). Multiplication of Ap, As and the reaction rate at theionization point gives the detected beam intensity. For assessment of observability in theexperiment, the intensity of the detected beam is evaluated in sHIBP. The detected beam
intensity must be large enough for observation, and is predicted to be larger than the noiselevel given by the detector.

FIG. 2. Variations along the heavy ion trajectory. (a) Attenuation rate of the primary andsecondary beams (ln Ai), (b) magnetic fields [T], (c) electrostatic potential [MeV], (d) electricfields [V/m], and (e) heavy ion velocity [106 m / s] are shown. The vertical dashed line
indicates the position where the ionization occurs from charge number +1 to +2.
The observation points depend on the injection conditions of the heavy ions, so the profile
scan can be made. The deviation of the potential values from the mean radial profile comesfrom 1, which depends on the observation point. 2-D profile scan is also possible withchange of the injected beam energy. Beam intensity includes effects of (1) damping in the
primary and secondary paths, and (2) ionization rate at the ionization point, so comparison
between the potential variation and beam intensity gives information on relative strength ofeffects (1) and (2).
Acknowledgements
This work was partially supported by JSPS KAKENHI JP17H06089, JP20K03905, and thecollaboration program of NIFS (NIFS22KIST029, NIFS23KIST049) and of RIAM of Kyushu
University. Some calculations were carried out using the JFRS-1 supercomputer system at
IFERC-CSC in QST and Plasma Simulator supercomputer system in NIFS.
References[1] LIN, L., et al., Phys. Plasmas 16 (2009) 012502.
[2] HOLLAND, C., et al., Phys. Plasmas 16 (2009) 052301.
[3] ROST, J. C., et al., Phys. Plasmas 17 (2010) 062506.
[4] KRUTKIN, O.L., et al., Nucl. Fusion 59 (2019) 096017[5] KUDASHEV, I., et al., Appl. Sci. 12 (2022) 9807.
[6] KASUYA, N., et al., Plasma Sci. Tech. 13 (2011) 326.
[7] KASUYA, N., et al., Nucl. Fusion 58 (2018) 106033.
[8] KASUYA, N., et al., Proc. BPSI meeting (2022) PA-5.
[9] YOSHIHARA, N., et al., Proc. BPSI meeting (2023) PA-5.
[10] IDO, T., et al., Rev. Sci. Instrum. 77 (2006) 10F523.
[11] FUJITA, K., et al., Nucl. Fusion 61 (2021) 086025.
[12] HIRSHMAN, S.P., et al., Comput. Phys. Commun. 43 (1986) 143.

Extraction of abrupt turbulence driven
particle transport using multi-ﬁeld SVD
T. Kodahara1, M. Sasaki1, Y. Kawachi2, Y. Jajima1, T. Kobayashi3, Y. Kosuga4,5,
G. Yatomi6, M. Nakata3,6,7, T. Yamada5,8, H. Arakawa9, A. Fujisawa4,5
Nihon Univ., 2KIT, 3NIFS, 4RIAM Kyushu Univ., 5RCPT Kyushu Univ.,
6
SOKENDAI, 7RIKEN iTHEMS, 8FAS Kyushu Univ., 9FMS Kyushu Univ.
1

1. Introduction

Transport in magnetized plasma is often accompanied by abrupt behaviors such as ELM andbreathers. Given the damage to the device wall caused by such transport, it is important tounderstand and predict the abrupt transport. However, since abrupt transport isspatiotemporally localized, Fourier analysis is generally unsuitable because it introduces manydegrees of freedom. Therefore, decomposition in fractional degrees of freedom by data-driven
scientiﬁc methods is easy to understand.
In this study, we propose the abrupt transport extraction method based on singular value
decomposition (SVD) [1, 2], which is one of the data-driven scientiﬁc methods. Density
ﬂuctuations and radial ﬂow ﬂuctuations evaluated by the linearly magnetized plasma
apparatus PANTA are used as targets. The radial ﬂow ﬂuctuation is evaluated from theelectrostatic potential by using the E×B drift process. Both ,luctuations are decomposed by
the same basis function using SVD, and the abrupt transport is evaluated from product of eachmode.
2. Multi-field SVD
The SVD is a diagonalization of a general matrix. Observables with spatio-temporal structure
can generally be described as φ(x, t); by using SVD for φ, it can be decomposed based on basis
functions as follows [1, 2]
φ = 𝑈𝛴𝑉 !

(1)

where Σ=diag(𝜎" , 𝜎# , …, 𝜎$ ), represented as a diagonal matrix (𝜎" >𝜎# >…>𝜎$ >0). The spatial
structure 𝑈$ and time evolution 𝑉$ corresponding to each 𝜎$ are also obtained. To extend
the SVD to multiple ,ields, we combine multiple physical quantities as follows.
𝑋 = [𝑓" (𝒙, 𝑡" ), ⋯ , 𝑓" (𝒙, 𝑡% ), 𝑓# (𝒙, 𝑡" ), ⋯ , 𝑓# (𝒙, 𝑡% ) ](2)

Using SVD for matrix 𝑋, the decomposition is based on common basis function; a multiple-,ield

SVD was used for the density ,luctuations 𝑁 and radial ,low ,luctuations 𝑉& obtained fromPANTA [3, 4]. By multiplying each singular value 𝜎$ 𝑈$ 𝑉$! , SVD modes 𝑁' and 𝑉&,) , which
represent the spatio-temporal structure at mode k, are obtained [5, 6].
The mode-by-mode product of both physical quantities yields the transport matrices Γ&,',) ,which represent the contribution to transport driven by each mode combination. Where α,βare the mode indices of the density ,luctuation 𝑁' in mode α and the radial ,low 𝑉&,) inmode β, respectively. By decomposing the transport matrix into a diagonal component Γ&,*+,- ,with equal and α and β other off-diagonal components Γ&,.// , we can decompose it into
stationary and abrupt components of transport. The spatio-temporal structure of both transportcomponents can also be discussed by summing over α and β as follows [7].
Γ&,*+,- = 5 𝑁' 𝑉&,'(3)

Γ&,.// = 5

(4)

'

'0)

𝑁' 𝑉&,)

3. Off-diagonal component transport
We discuss the off-diagonal components of the transport matrix. For the spatio-temporal
structure of the off-diagonal component transport, a two-dimensional Fourier mode
decomposition reveals a broad spectrum in both the azimuthal and frequency directions. Thismeans that short-lived structures appear in time and space. In addition, the spectrum extends inthe positive wavenumber direction, indicating that it propagates in the positive direction. This
corresponds to the electron-antimagnetic direction, which is consistent with the direction ofoscillation propagation. In addition, the time evolution at a certain azimuthal point showsintermittent large amplitude transport. Compared to the time average of the diagonalcomponent in the transport matrix, the transport is about 100 times more frequent.
Next, the behavior of the oscillations when abrupt transport occurs in the off-diagonalcomponent is investigated. From the transport matrix, we extract the modes that contribute to
the abrupt transport and compare the temporal and spatial behavior of the abrupt-driven modeand the off-diagonal component transport in both oscillations. Figure 1 (left panel) shows thetime evolution of the off-diagonal component transport and both ,luctuations. In the off-diagonalcomponent, when abrupt transport occurs, the amplitude of both ,luctuations is large and inphase with each other. The duration of the abrupt transport is about 1/10 of the ,luctuation
period and occurs instantaneously. The spatial structure of the abrupt transport when it occurs
is also compared with that of the ,luctuation. Figure 1 (right panel) shows the spatial structureof the off-diagonal component and both ,luctuations. In the spatial structure, it can be seen thatthe amplitudes both ,luctuations are large and in phase with each other when abrupt transport
occurs. In addition, the peak width of the abrupt transport is about 1/5 of the wavelength of theturbulence and is spatially localized. Thus, the decomposition into nonlinear modes using SVD

%!
&",$

#",$

!!

!!

%!
&",$

#",$

"" !"""" !""

allows us to identify such modes that drive the abrupt transport.

Figure 1 time evolution (left panel) and spatial structure (right panel) of off-diagonalcomponent transport and both ,luctuations. The red line represents the density and the blueline the radial ,low. The amplitudes of both ,luctuations are large in time and space, and abrupt
transport occurs at the moment of phase coincidence.
Reference[1] M. Sasaki, et al., Plasma Phys. Control. Fusion, 63.2, 025004 (2020).
[2] Y. Asahi, et al., Physics of Plasmas, 28.1, (2021)[3] S. Inagaki, et al., Scienti,ic reports, 6.1, 22189 (2016).
[4] T. Kobayashi, et al., Plasma Fusion Res., 12, 1401019 (2017).
[5] G. Yatomi, et al., Plasma Phys. Control. Fusion, 65.9, 095014 (2023).
[6] T. kodahara, et al., Plasma Fusion Res., 18, 1202036 (2023)[7] F. Kin, et al., Physics of Plasmas, 26.4 (2019)

Optimizations of NB deposition calculations in the integrated code
*Kazutoshi Kani, Takaaki Fujita, Atsushi Okamoto, Sota Takemoto, Ren Kadomatsu
Nagoya Univ.

1

Introduction
The calculation speed of integrated codes needs to be more accelerated for the application to real time control

and to surveying parameters. As for TOTAL code utilized in our group, calculation for NBCD (Neutral Beam
Current Drive) takes long time. Approximately 60% of the calculation time is used for the NBCD calculation.
Neutral particles injected with NB into the plasma become fast ions through interaction with electrons andions in the plasma, orbiting the torus, resulting in a plasma current. The conceptual diagram is shown in
Fig.1. In this study, we optimized NBCD calculation through two approaches for faster code executing. First,
investigation of optimal number of test particles. Second, modify the code to generate multiple fast ions along
a single trajectory of beam atoms.

Fig.1 Conceptual diagram of NBCD

2

NBCD calculation methods
The conceptual diagram of determinating fast ion generation poiont is shown in Fig.2. The procedure ofcalculating NB current are the followings.
1. Finding fast ion generation point.
i. Randomly generate a test particle on the cross-section of the ion source.(Cross-section of the ion
source is divided into 200 sections in both r and θ derections. )ii. Calculating the non-ionization probability of a test particle along the beamline.
iii. The point of fast-ion generation is set when the non-ionization probability falls below the ionizationthreshold (red line in Fig.3) determined by a random number in advance.
iv. Repeating i～iii for 104 times.
2. Determine the fast ion sourse.
3. Solving Fokker-Plank equation to obtain NB current drive density distribution.
Since the conditions for calculating the non-ionization probability, such as plasma temperature and density,
depends on time and ion source cross-section location, the non-ionization probability must be recalculated eachtest particle.

1

Fig.2 Conceptual diagram of fast
ion generation point determination
along the beamline.

3

Fig.3 Conceptual diagram of finding the fast-ion generating point.
The point where the blue line andthe red line intersect is the point of
fast-ion generation.

Fig.4 NB driven current density
distribution at t = 100 s.

Optimization 1
The first optimization for speeding up the code is investigation of the optimal number of test particles. The

way to investigate is as follows. First, run code with fewer test particles generated. Second, calculate the error
rate in the NB-driven current density at each radial position. The error rate is defined as follow.

(error-rate[%]) =

jN B

104

−jN BnjN B

104

× 100 , where jN Bn (104 ) denotes NB current density in test particles

n(10 ). Then average computation time per step during the period with small change in plasma parame4

ters. Finaly determine optimal number of test particles based on error rate and reduced computation time.
The results are shown in Fig.5. Fig.5 shows that there is approximate linear relationship between the numberof particles generated and average computation time per step. We can reduce NB computation time by 65% by
reducing the number of test particles to 1/10. And even if the number of test particles is reduced to 1/10, theerror rate is about 4%.

Fig.5 Relation between test particles and one step calculation time.

2

4

Optimization 2
The second optimization for speeding up the code is modify the code to generate multiple fast ions along

a single trajectory of beam atoms. In a conventional method, firing one test particle from one point on the
cross-section of the ion source to find one fast ion generation point at a step. We modified the code to set n test
particles emitted from a single cross section point and to continue the calculation until all test particles havebeen ionized. (or, until the test particles penetrate the plasma.) The conceptual diagram of this new method
is shown in Fig.6.

Fig.6 Conceptual diagram of optimization 2 (n = 3).
Modified the code to allow n test particles to be fired fromone point and n fast ion generation points to be obtained
at a step.

Fig.7 Relationship between thenumber of test particles per one test
particle and the computation time
per step for n = 1 and n = 5.

Figure 7 shows the relationship between the number of test particles per one test particle and the computation
time per step for n = 1 and n = 5. We confirmed that applying this optimization could reduce calculation time.
However, the appropriate value of n, or the possibility that this logic is inappropriate (since the computation
time is longer in the new method with n = 1 than in the conventional method), is an issue for future reserch.

5

Summary and Prospect
In summary, we attempted to optimize the NB calculation using two approaches. First, we invested the

optimal numbers of test particles. The results shows that reducing the number of test particles to 1/10 reduces
NB computation time by 65%. Second, modifying the code to generate multiple fast ions along a single trajectoryof beam atoms.
For the future prospects, we will determine the appropriate number of test particles emitted from a single
cross section point. Also, parallelizing of non-sequential computation points is being considered.

3

Simulation of PCI experimental diagnostic for density fluctuations
in tokamak plasmas
Kosuke KUWAMIYA1, Naohiro KASUYA1,2, Shota MOCHINAGA1
IGSES, Kyushu Univ.1, RIAM, Kyushu Univ.2

1.

Research background and objectives
Improving the confinement performance of plasma is crucial for the realization of nuclear fusion power

generation. Since turbulent transport plays a significant role in the plasma confinement [1], understanding the
turbulent transport is important. Phase Contrast Imaging (PCI) is one of the methods for measuring turbulence inmagnetically-confined plasmas. For PCI, phase change of an injected laser beam is measured due to density
fluctuations [2,3]. As PCI does not depend on specific plasma conditions and allows measurement across the entire
plasma region, it has been adopted in various tokamaks and helical devices. Numerical simulations have also been
conducted for comparison with experiments [2,4].
In this study, we perform PCI measurement simulations for experimental device JT-60SA. We simulate PCI
measurements by combining magnetic field equilibrium calculations in the JT-60SA tokamak with fluctuation fields
generated by a reduced MHD code as test data. Evaluation of the line integration effect in the direction of the lineof sight for density fluctuations gives understanding of the signals obtained from PCI measurements.

2.

Turbulence measurement simulation
We perform simulations of PCI measurements on fluctuation field data in a three-dimensional tokamak magnetic

configuration. This section describes an overview of that process.
2.1 Plasma Equilibrium
For the simulation of measurements, the geometry of the device and the shape of the plasma are necessary. A
tokamak plasma shape is calculated by using the integrated code TASK/EQU module. The EQU module is capableof performing equilibrium calculations for free-boundary plasmas, including the external coils currents [5].
2.2 Turbulence field
The turbulence to be measured is obtained by turbulence simulations. Here, we use pressure-driven ballooning
type fluctuations calculated by using the reduced MHD code R5F [6]. Time evolution of five fields—potential,
vector potential, parallel flow velocity, density, and electron temperature—are solved spatially in three dimensions
using the vorticity equation, Ohm's law, the parallel ion motion equation, the continuity equation, and the electron
temperature evolution equation. We use these fields as test data, although effects of the equilibrium shape are notincluded in the turbulence simulation.

2.3 PCI measurement simulation
The PCI measurement gives the wavenumber spectrum of density fluctuations [3]. The incident laser is scatteredby the density fluctuations, and the interference of this scattered wave with the reference wave provides a twodimensional pattern related to the wavenumber k of the fluctuation. Since the observation signal is integrated one
along the line of sight, it becomes necessary to reconstruct local density fluctuation information from this signal.
Here images of the density fluctuations integrated along the measurement line of sight are provided as themeasurement results.

3.

Results of data analysis

3.1 Equilibrium calculation
Figure 1 shows the magnetic field configuration obtained by using TASK/EQU. The coil currents used for the
calculation are also shown in Table 1. In this case, the co-direction to the plasma current gives a positive value. The
used parameters are major radius R = 3.0m, minor radius a = 1.2m, plasma volume Vp = 133m³, and toroidal
magnetic field at the plasma center BT = 2.25T, as for JT-60SA. As a result, an equilibrium magnetic field
configuration with triangularity δ = 0.38 and elongation κ = 1.85 is obtained.
Table. 1 coil current

Fig. 1 Equilibrium magnetic field configuration

3.2 PCI measurement simulation
On the magnetic field equilibrium obtained by EQU, turbulence simulation data obtained by R5F is placed as
the three-dimensional distribution of fluctuations (Fig. 2 (c)). In JT-60SA, the laser is injected tangentially to themagnetic field direction (Tangential Phase Contrast Imaging, TPCI) [2].
Using the TPCI beam path, as shown in Fig. 2, we calculate the integrated density fluctuation image. In JT-60SA,
there is difference in the vertical direction of the injection and detection positions [2]. Therefore, in this simulation,
we carry out the simulation along the path not horizontal. Figure 3 shows the snapshots of local (at the position B3
in Fig.2) and integrated patterns of density fluctuations, and the calculated two-dimensional wavenumber spectrum
from them. The axis x and y in Fig. 3 correspond to the direction perpendicular and horizontal to the line of sight,
respectively.
Focusing on the two-dimensional wavenumber spectrum of the integrated one, peaks are observed at (kxρi, kyρi)= (± 0.54, ± 0.18). Converting kxρi = ± 0.54 to wavelength gives λx ≈ 0.035 m, which matches the wavelength of thefluctuation in the x-direction of the local data. This indicates that fluctuations of this scale appear in the spectrum.
On the other hand, converting kyρi = ± 0.18 to wavelength gives λy ≈ 0.10 m, which corresponds to the beam width(0.10 m). Actually, the wavelength of the dominant component in the y-direction is 0.20 m, so is not possible to beresolved in this measurement.
(a)

(b)

(c)

Fig. 2 TPCI measurement line of sight (orange lines). The views from (a) above and (b) the side are shown,with B1 representing the laser incidence position and B2 the detection position. The fluctuations andthe line of sight on the poloidal cross-section are also shown in (c).

(a) density fluctuation

wavenumber spectrum

local

local

(b) density fluctuation

wavenumber spectrum

integrated

integrated

Fig. 3 Density fluctuations. Local and integral images and their wavenumber spectra are shown.

4.

Conclusion
In this study, we constructed a PCI measurement simulation routine and performed calculations using themagnetic field configuration of JT-60SA tokamak plasma. We applied the TPCI diagnostic arrangement as in JT60SA, and presented the local and line-of-sight integrated density fluctuations, as well as their two-dimensionalwavenumber spectrum.
In the future, we will develop a routine to construct local information from the integrated density fluctuation
data. We can utilize the fact that the direction of the magnetic field lines changes in relation to the measurement lineof sight. Moreover, by analyzing microscopic turbulence simulation data considering the magnetic field
configuration, we aim to deepen our understanding of turbulent transport caused by plasma instabilities.

References[1] P. H. Diamond, et al., Plasma Phys. Control. Fusion 47, R35 (2005).
[2] S. Coda, et al, Nucl. Fusion 61, 106022 (2021).
[3] K. Tanaka, et al., Rev. Sci. Instrum. 79, 10E702 (2008).
[4] N. Kasuya, et al., Nucl. Fusion 58, 106033 (2018).
[5] 松岡啓介，九大応力研所報 141 51 (2011).
[6] M. Yagi, et al., 26th IAEA FEC, TH/P3-21 (2016).

Evaluation of nonlinear mechanism to induce inward turbulent
particle fluxes in tokamaks
R. Miyamoto1, N. Kasuya1,2, M. Yagi3
1

IGSES Kyushu Univ., 2RIAM Kyushu Univ., 3QST.

1. Introduction
In a fusion reactor, identifying an effective method for particle supply to core plasma is an
important issue to maintain high density plasma. The pellet injection method is one of the ways tosupply fuels to the core. With pellet injection, peaked density profiles are transiently formed andassociated fluctuations have been observed [1]. In such regions, it is known that a particle flux
toward the center is generated by fluctuations [2]. Therefore, a global turbulence simulation was
conducted to evaluate the nonlinear mechanism of density fluctuations, which drive particle fluxes[3].
2. Simulation model
The reduced MHD code R5F [4] is used for the simulations. The code solves the time evolutionof 3-D fields (electrostatic potential 𝜙, vector potential in the magnetic field direction 𝐴, velocity inthe magnetic field direction 𝑉, density 𝑁, and electron temperature 𝑇𝑒 ) using the following
vorticity equation, ohm's law, equation of ion motion in the magnetic field direction, continuity
equation, and evolution equation for electron temperature;
𝑑 2
𝑞 1𝜕
𝑞 𝜕𝜙
𝑞 𝑚𝑒 1 𝜕𝑞 𝜕(𝜙 − 𝛿𝑝𝑒 )∇ 𝜙 = −∇∥ ∇2⊥ 𝐴 + 𝜇𝑖⊥ ∇4⊥ 𝜙 − [2𝑟𝑐𝑜𝑠𝜃, 𝑝] − 𝜇𝑖𝑟 {𝑣 +}−
𝜇𝑟 {𝑣 + 𝛿𝛻⊥2 𝐴 +},𝑑𝑡 ⊥𝜀 𝑟 𝜕𝑟𝜀 𝜕𝑟𝜀 𝑚𝑖 𝑒 𝑟 𝜕𝑟
𝜀
𝜕𝑟

(1)

𝑑
𝑚𝑒 2
𝑚𝑒 2 2𝑚𝑒𝑞 𝜕(𝜙 − 𝛿𝑝𝑒 )(𝐴 − 𝛿 2∇ 𝐴) = − ∇∥ (𝜙 − 𝛿𝑝𝑒 ) + 𝜂∥ ∇2⊥ 𝐴 − 4𝜇𝑒⊥ 𝛿 2
𝛻 𝛻 𝐴 + 𝛼𝑇 𝛿𝛻∥ 𝑇𝑒 + 𝛿
𝜇 {𝑣 ++ 𝛿𝛻⊥2 𝐴} , (2)𝑑𝑡𝑚𝑖 ⊥𝑚𝑖 ⊥ ⊥
𝑚𝑖 𝑒
𝜀
𝜕𝑟
𝑑𝑉
𝑞 𝜕𝜙
𝑚𝑒
𝑞 𝜕(𝜙 − 𝛿𝑝𝑒 )= −∇∥ 𝑝 + 4𝜇𝑖⊥ ∇2⊥ 𝑉 − 𝜇𝑖 (𝑉 +)−𝜇 {𝑉 ++ 𝛿𝛻⊥2 𝐴} ,
𝑑𝑡
𝜀 𝜕𝑟
𝑚𝑖 𝑒
𝜀
𝜕𝑟

(3)

𝑑𝑁
𝑑𝑝
+𝛽
= −𝛽∇∥ (𝑣 + 𝛿∇2⊥ 𝐴) + 𝜂⊥ 𝛽∇2⊥ 𝑝 + 𝛽[2𝑟𝑐𝑜𝑠𝜃, 𝜙 − 𝛿𝑝𝑒 ],
𝑑𝑡
𝑑𝑡

(4)

3 𝑑𝑇𝑒 𝛽𝑒 𝑑𝑁
5−= −𝛼𝑇 𝛿𝛽𝑒 ∇∥ ∇2⊥ 𝐴 + 𝜒𝑒⊥ ∇2⊥ 𝑇𝑒 + 𝜖 2 𝜒𝑒∥ 𝛻∥2 𝑇𝑒 − 𝛿𝛽𝑒 [2𝑟𝑐𝑜𝑠𝜃, 𝑇𝑒 ],
2 𝑑𝑡
𝛽 𝑑𝑡
2

(5)

(0)

where 𝑑𝑡 ∶= 𝜕𝑡 + [𝜙, ], ∇∥ ∶= ∇∥ − [𝐴, ]. The set of equations is normalized by using the minor
radius 𝑎 for the length and by the poloidal Alfvén time 𝑎/(𝜀𝑣𝑎 ) for the time, where 𝜀 is theinverse aspect ratio and 𝑣𝑎 is the Alfvén velocity. The calculation parameters used for thesimulation are as follows; major radius 𝑅=1.5 m, minor radius 𝑟=0.5 m, beta value 𝛽 = 𝛽𝑒 =5 × 10−3, inverse aspect ratio 𝜀 = 0.33, ion skin length 𝛿 = 10−2, mass ratio 𝑚𝑒 ⁄ 𝑚𝑖 =

1/1836, resistivities 𝜂⊥ = 10−5, 𝜂∥ = 10−5 , thermal diffusion coefficients 𝜒𝑒⊥ = 10−5 , 𝜒𝑒∥ =10−5, neoclassical viscosity coefficients 𝜇𝑒⊥ = 10−5 , 𝜇𝑖⊥ = 10−5, and thermal coefficient 𝛼 𝑇 =0.71. In this study, the calculations are carried out with an initial density profile peaked at radial
position 𝑟𝑠 = 0.8.
3. Analysis results
The initial density profile is set to be 𝑁𝐸𝑄 = 𝑁0 (1 − 𝑟 8 )8 + 𝑆𝑎𝑚𝑝 exp(− 𝜉 2 ⁄2𝛥2) with 𝜉 =𝑟 − 𝑟𝑠 , 𝑟𝑠 = 0.8, 𝛥2 = 0.003, and the cases with 𝑆𝑎𝑚𝑝 = 1.0 and 1.5 are compared [3]. Figure 1
shows the radial profiles of the particle flux for 𝑆𝑎𝑚𝑝 = 1.0 and 1.5. In both cases, inward particle
fluxes were observed around 𝑟 = 0.75, where the safety factor is q ≅ 2.6*. Figure 2 shows the time
evolution of the particle flux at 𝑟 = 0.75. The inward particle flux is about 10 times larger and issustained for a longer duration for 𝑆𝑎𝑚𝑝 = 1.5 than 𝑆𝑎𝑚𝑝 = 1.0.

Fig.1 Radial profiles of the particle flux in theFig.2 Time evolution of the particle flux

nonlinear states for 𝑆𝑎𝑚𝑝 = 1.0 and 1.5.

for 𝑆𝑎𝑚𝑝 = 1.5 and 1.0 at 𝑟 = 0.75

Positive values correspond to the radiallyoutward direction.
The turbulent particle flux 𝛤𝑡𝑢𝑟𝑏 is evaluated by the following equation from the nonlinear
convective derivative term in the continuity equation.
𝛤𝑡𝑢𝑟𝑏 = ∑𝑚,𝑛{(1 + 𝛽)𝑁𝑚,𝑛𝜕𝜙−𝑚,−𝑛𝜕𝑣−𝑚,−𝑛𝜕∇2⊥ 𝐴−𝑚,−𝑛+ 𝛽𝐴𝑚,𝑛+ 𝛽𝛿𝐴𝑚,𝑛
}
𝑟𝜕𝜃
𝑟𝜕𝜃
𝑟𝜕𝜃

(6)

Therefore, (𝑚, 𝑛) ≠ (0,0) modes must be sustained for particle flux generation, where 𝑚 is thepoloidal mode number and 𝑛 is the toroidal mode number. In this study, we focused on the density
fluctuation of the first term in Eq. (6) because its contribution to the particle flux is relatively larger
than other two terms. Figure 3 shows decomposition of the particle flux into contribution of eachdensity fluctuation component. There are relatively low wavenumber modes, such as (𝑚, 𝑛) =(−3,1), (−5,2), and (−8,3), which have larger contribution to the particle flux in the case 𝑆𝑎𝑚𝑝 =
1.5.

(𝑎)

(𝑏)

(−5,2)(−8,3)(−3,1)

Fig.3 Decomposition of the particle flux into contribution of each density
fluctuation components for cases (𝑎)𝑆𝑎𝑚𝑝 = 1.0 and (𝑏) 𝑆𝑎𝑚𝑝 = 1.5
Figure 4 shows the nonlinear contributions to the (−3,1) density fluctuation (𝑁−3,1) for 𝑆𝑎𝑚𝑝 =1.0 and 1.5. Here, the averaged values for 𝑡 = 113~137 and 𝑡 = 72~190 are plotted for the
𝑆𝑎𝑚𝑝 = 1.0 and 𝑆𝑎𝑚𝑝 = 1.5 cases, respectively, corresponding to the time period when the inward
fluxes are generated. Figure 4(a) shows that the nonlinear contribution of higher wavenumber
components around |𝑚| ~ 20 is significant in the 𝑆𝑎𝑚𝑝 = 1.0 case. On the other hand, Figure 4(b)shows that for 𝑆𝑎𝑚𝑝 = 1.5, there are various frequency components including lower wavenumber
components that have significant influences, and their magnitudes are more than 10 times larger thanthose for 𝑆𝑎𝑚𝑝 = 1.0.
(𝑎)

(𝑏)

Fig.4 Nonlinear contributions to the (-3,1) density fluctuation for (𝑎)𝑆𝑎𝑚𝑝 = 1.0
and (𝑏)𝑆𝑎𝑚𝑝 = 1.5. The averaged values are plotted for the time period wheninward particle fluxes are generated.

Next, the nonlinear contributions for the 𝑆𝑎𝑚𝑝 = 1.5 case are evaluated. Figure 5 shows that
nonlinear contributions to (𝑎) (−3,1), (𝑏) (−5,2), (𝑐) (−8,3), and (𝑑) (−11,4) modes, which
are rather lower wavenumber modes contributing to drive particle fluxes. Those are the averagecomponent for 𝑡 = 72~190. In this time period, the larger particle flux is sustained for 𝑆𝑎𝑚𝑝 =1.5, compared with 𝑆𝑎𝑚𝑝 = 1.0. From Fig. 5, the following three-wave couplings are found to havesignificant contributions;(−3,1) ← (−8,3) + (5, −2)(−5,2) ← (−11,4) + (6, −2)(−8,3) ← (−11,4) + (3, −1)(−11,4)← (−14,5) + (3, −1){(−3,1)(𝑎)(5, −2)(−5,2)(𝑏)(−11,4)(−8,3)(𝑐)(6, −2)(3, −1)(𝑑)(3, −1)

Fig.5 Nonlinear contributions to (𝑎) (−3,1) , (𝑏) (−5,2) , (𝑐) (−8,3) , and(𝑑) (−11,4) modes for 𝑆𝑎𝑚𝑝 = 1.5.
Therefore, in the case of 𝑆𝑎𝑚𝑝 = 1.5, the fluctuation is sustained by the low wavenumber modes
having a three-wave coupling relationship as in Fig.6, which induces a relatively larger inward
particle flux for 𝑆𝑎𝑚𝑝 = 1.5 than for 𝑆𝑎𝑚𝑝 = 1. 0.

( 3,1)( 5,2)( 8,3)( 11,4)

Fig.6 Three-wave coupling relationship for 𝑆amp = 1.5.
4. Conclusion
The R5F code was used to analyze the nonlinearly driven particle flux by electron drift mode.
Nonlinear couplings are decomposed into contributions of each Fourier component. With thestronger density gradient case, when an stronger inward particle flux is driven, it is found that
relatively lower wavenumber modes as (𝑚, 𝑛) = (−3, 1), (−5, 2), (−8, 3) have large contributions
and the fluctuations are sustained by their mode couplings.
Reference[1] S. Ohshima, et al, Sci. Reps., 12 (2022) 14204.
[2] M. Yagi, et al, Proc. 16th BPSI Meeting (2018) 2-1.

[3] H. Todoroki, et al, Plasma Fusion Res. 18 (2023) 1203052.
[4] M. Yagi, et al, 26th IAEA FEC (2016) TH/P3-21.

Analysis of parallel flow shear driven mode including
effect of finite ion temperature gradient
Itsuki OYAMA1 and Yusuke KOSUGA2
1

IGSES, Kyushu University

2

1

RIAM, Kyushu University

Introduction

A good confinement requires transport barriers, for example external transport barrier (ETB)
and internal transport barrier (ITB). As the experimental facts, it has been identified that
there is such transport barriers in the experimental devices[1]. When transport barrier is built,
plasma profiles get several kinds of gradients. Several mechanisms for losing transport barrier,especially ITB, are still missing. As the first step of understanding the cause to lose ITB, it is ofimportance to pay attention to a turbulent transport. This is because turbulent transport canaffect the stability of ITB through some effect to the profile.
To shed light on it, we investigated how instability growth in the situation where three
gradients (density, ion temperature and parallel velocity) coexist. It must be noted that parallel
velocity gradient means the gradient of velocity parallel to the toroidal magnetic field lines.

2

Model

In this study, 3 gradients (density, ion temperature and parallel flow) are assumed. Especially
we paid attention to the coexistence relation between ITG mode and PVG (Parallel Velocity
Gradient driven) mode. To elucidate it, we use the set of fluid equations shown below :



∂ eϕ̃
∇n0
c d 2
pi



+ vE ·+ ∇∥ ṽ∥ −
∇
ϕ̃ +
=0


∂t Te
n0
Bωci dt ⊥en0∇∥ p̃ideṽ∥ + vE · ∇V∥ = −
−
∇ ϕ̃

dt
mi n0 mi ∥






5
3 ∂
5

(2)


+ vE · ∇ (P0 + p̃) + P0 ∇∥ ṽ + P0 ∇ · v⊥ + ∇ · q (2) = 0
2 ∂t
2
2

(1)

where ...
˜ is a fluctuation. Subscripts, i or e, means ion and electron. e is the elementary electriccharge. ϕ is the electrostatic potential. p is the pressure composed of density n and temperature1T . B is the magnetic field. ωci is the ion cyclotron frequency. First equation comes from the lawof density conservation. Second one comes from the momentum equation in a parallel direction.
Third equation comes from the law of energy conservation. Last two terms in third equation are
already reported in Pogtsue et al.,[2]. It is to be noted that we assume the adiabatic response
between electron density perturbation ñ/n0 and electrostatic potential eϕ̃/Te0 .

3

Analysis

To elucidate which instability can be destabilized most, we derived the dispersion relation ofthis set of equations.
!"
#




2 2
10 2 210 2 2ω∗i
5 k∥ ci10
2 2
2 2 ω∗pi
1 + ρi k⊥ −1+ τ +
ρi k⊥ − 1 + ρi k⊥ τ
+ ρi k⊥ τ
3
3 ω2
3
3
ω
ω






k∥2 c2i ω∗pi
1 5 2 25 2 2
10 2 2 ky ρi′ρ k +− 1 + ρi k⊥⟨V ⟩= 2 τ 1 + ρi k⊥ 1 + 1 +ω3
τ 3 i ⊥
ω3k∥ ci ∥

(2)

where τ is the temperature ratio Te /Ti . ω∗pi = (1 + ηi )ω∗i . ηi is the ratio between density
gradient and ion temperature gradient d ln Ti /d ln n and ω∗i is the ion drift frequency.
As the first step of analysis, we supposed the long wavelength. This assumption leads to the
reduction of Eq. (2) by using ρi k⊥ −→ 0. Then, Eq. (2) can get formed as :2 2 

ω∗i k∥ ci
− 2
1−
ω
ω




 
ky ρi
5
2 ω∗i
′
=0
+ 1−
⟨V ⟩ + ηi −3k∥ ci ∥
3 ω

(3)where τ is assumed as 1 for this time. The situation τ ∼ 1 corresponds to the situation where
ion and electron are in almost same temperature. To get the growth rate of instabilities, we
calculated Eq. (3) numerically :Figure 1: ITG mode dominated region.

Figure 2: PVG mode dominated region

2

Because Eq. (3) is obviously cubic, there are three solutions for this. However one of three has
negative growth rate (γ ≤ 0), so it is omitted. In these graphs, the vertical axis corresponds toηi = d ln Ti /d ln n and horizontal axis corresponds to the strength of parallel velocity gradient⟨V∥ ⟩′ /(cs /Ln ). As it goes up in the vertical axis, ion temperature gradient is developed. As forthe horizontal axis, parallel velocity gradient is developed most at both end.
In Fig. 1, as ion temperature gradient is surged, it is found that the value of growth rate
goes higher. Then, we can recognize this as ITG mode. Paying attention to the lowest curve(γ = 0), the threshold of ITG mode goes higher, as velocity gradient is developed. This is why
parallel flow shear stabilizes ITG mode. When parallel flow shear becomes strong enough to get
ITG mode suppressed, we can see that situation in Fig. 2.
In Fig. 2, as parallel velocity gradient is developed, it is found that the value of growth rate
goes higher. Then, we can recognize this as PVG mode. When ηi gets developed, the growthrate goes lower. This is why finite ion temperature gradient also stabilizes PVG mode.
Following the long wavelength assumption, we investigated the short wavelength scale. Although analysis of the long wavelength modes is important to understand general trend of
instabilities, they are strongly affected by transport barriers. This is why analysis of the shortwavelength is more important on condition that transport barrier exists. To calculate the growthrate, we utilized Eq. (2).

Figure 3: Short wavelength ITG mode domi- Figure 4: Short wavelength PVG mode dominated region.

nated region

Here τ is still equal to 1. Compared Fig. 1 with Fig. 3, we can find the stabilization effects byparallel flow shear get strengthened. In addition, the ITG-damped region gets broader. In thesame way as long wavelength, PVG mode gets destabilized in ITG-damped region. However thiscontour says the stabilization effects by ηi is reduced in comparison with the long wavelength
case. Hence it can be said that the importance of PVG mode gets increased relatively.
3

Meanwhile it is also of note to consider the case where τ is finite. That situation is not
so strange, because heating plasma techniques, e.g. Electron Cyclotron Resonance Heating(ECRH), can help τ to get finite. Furthermore, cold ion assumption has already been recognized,
which ion temperature at edge plasma can be regarded as 0. In this case, τ is actually infinite.
We can get the contour map with finite τ by using Eq. (2) below.

Figure 5: Short wavelength ITG mode domi- Figure 6: Short wavelength PVG mode dominated region. τ = 10.

nated region. τ = 10.

Compared Fig. 1 with Fig. 5, one easily sees that ITG mode gets nearly damped. Obviously theregion, where ITG mode excites, gets shrunk, and the stabilization effect by parallel flow shear
gets stronger in this region. Moreover the value of growth rate gets smaller even at most. Due tothis, it is perhaps conceivable that ITG mode can not play a role as a source to drive transport
in this case. In contrast, PVG mode get its region spread and significantly destabilized in that
region. In that region, the stabilization effects by ηi is mostly gone. In consequence, it can beconsidered that PVG mode plays an important role as a driving source of transport in this kindof case.
This work is partly supported by the Grants-in-Aid for Scientific Research of JSPS of Japan(JP21H01066), the joint research project in RIAM, Kyushu University.

References[1] K. Ida et al., Plasma Phys. Control. Fusion (2018) 60 (3), 033001[2] I. O. Pogtse et al., J. Plasma Phys. 60 (1), 133-149 (1998)

4