
Citation: | Xiaotao XIAO, Shaojie WANG, Lei YE, Zongliang DAI, Chengkang PAN, Qilong REN. BITS: an efficient transport solver based on a collocation method with B-spline basis[J]. Plasma Science and Technology, 2023, 25(2): 025103. DOI: 10.1088/2058-6272/ac909f |
A B-spline Interpolation Transport Solver (BITS) based on a collocation method is developed. It solves transport equations as a generalized interpolation problem, taking the first-order accuracy in time and the second-order accuracy in space along with a predictor–corrector or under-relaxation iteration method. Numerical tests show that BITS can solve one-dimensional transport equations for tokamak plasma more accurately without additional computation cost, compared to the finite difference method transport solver which is widely used in existing tokamak transport codes.
The complexity of transport process in tokamak plasma is generally acknowledged due to various kinds of sources and sinks, instabilities, plasma–surface interactions and so on [1–3]. Numerical simulation is an essential tool to understand the transport physics of tokamak plasma. Since the magnetic field in a tokamak has a nested structure, the magnetohydrodynamic equation of particle balance, electron and ion energy conservation, Faraday's law, toroidal momentum can be simplified to a set of one-dimensional (1-D) flux surface averaged transport equations [4, 5]. A lot of transport codes have been developed by solving these transport equations since around the 1980s, such as ONETWO [6], TRANSP [7, 8], ASTRA [9], and ETS [10]. These codes are widely used to interpret experiments at modern tokamaks, such as EAST [11] and HL-2A [12], and help to predict transport of future tokamak devices, such as CFETR [13]. In order to simulate all phases of a tokamak scenario, the transport codes are integrated with other tokamak-physics codes. A lot of integrated suites of codes have been developed, such as CRONOS [14], JINTRAC [15], OMFIT [16], IMAS [17] and TRIASSIC [18].
Due to the huge computational resources required by these simulations, it is meaningful to develop an efficient transport solver. Although transport solvers of these codes are firstly and usually developed based on the finite difference method (FDM), different numerical methods continue to be developed to make the transport solver more efficient. For example, the finite element method using a Hermite basis function is applied to solve the current diffusion equation in tokamaks [19]. The finite volume method is used in a two-dimensional (2-D) transport solver as it guarantees conservation of the fluxes in both parallel and perpendicular directions [20]. Since widely used theory-based transport models such as GLF23 [21], MMM95 [22] and TGLF [23] have a strong nonlinear dependency on the temperature and density gradients, it is very important to obtain the gradients of physics variables accurately and smoothly. A high-order interpolated differential operation scheme [24] is employed to solve not only the physics variable but also its gradient. It provides an accurate transport solution using a small number of grid points with robust nonlinear convergence. However, this scheme is still under the framework of FDM; the result is improved at the cost of much more computation cost (the gradient of the variable is taken as additional independent variable and a higher-order accuracy formula is used). In this work, a B-spline Interpolation Transport Solver (BITS) is developed to solve the transport equations of tokamak plasma. The computation costs of solving differential equations by FDM and a collocation method are at the same level, while smoother results can be obtained by the collocation method. This helps to reduce the error of numerical results greatly.
In section 2 the transport equations, numerical scheme of the collocation method and FDM, and iteration method are introduced. The results of a few different cases using FDM and the B-spline collocation method are compared in section 3. A conclusion is given in the final section.
The 1-D energy diffusion equation for magnetically confined toroidal plasma in an axisymmetric tokamak can be simplified in cylindrical coordinates as
32n∂T∂t-1ρ∂∂ρ[ρ(nχ∂T∂ρ)]=S | (1) |
Here, n is the particle density, T is the temperature of the electron or ion species, ρ is the normalized minor radius proportional to the square root of the toroidal flux, χ is the effective thermal diffusivity, and S is the source term that includes ohmic heating, ionization, radiation loss, electron–ion heat exchange, and external auxiliary heating such as neutral beam injection and radiofrequency heating. For a given source S(ρ, t), boundary condition
The differential equation can be solved as a generalized interpolation problem. In a classical interpolation problem, the constraint equation is that the combination of coefficients αj and basis functions Bj(x) is equal to the given value yi at points xi, i.e. ∑jαjBj(xi) = yi. When solving the differential equation, the constraint is that the combination of coefficients αj, basis functions Bj(x) and their derivatives
Based on collocation, T is approximated by a linear combination of basis function Bi, k(ρ). Here the second-order accuracy is used, i.e. T(ρ) = ΣiαiBi, 3(ρ). The stable evaluation of B(ρ) and its derivatives are given in the appendix. It can be simply understood as a generalized piecewise polynomial function. It is nonzero only in k successive intervals, and it is actually a polynomial function on each interval. Once grid points are given, and a proper knot sequence t1, ···, tN+k is generated, then Bi, k(ρ) is totally determined. In this work, the knot sequence is generated as follows:
(2) |
Then temperature T at grid points
(3) |
in which MBρ is a tridiagonal matrix, with only two additional nonzero elements m1, 3, mN, N-2. Here element MBρ(i, j) is written as mi, j = Bj(ρi).
(4) |
As temperature T varies with time, so does the coefficient. We have
(5) |
(6) |
in which
(7) |
(8) |
in which
(9) |
Here MBρ,
For interior points i = 2, ···, N - 1, the equation is discretized as
(10) |
For the inner boundary condition
(11) |
For the outer boundary condition
(12) |
It is easy to construct a general boundary condition in a similar way. For example,
(13) |
Finally, equations (10)–(12) constitute a matrix equation,
(14) |
in which MCM is a tridiagonal matrix, with two additional nonzero elements (the third one in the first row and the third from last one in the last row);
It is very easy to implement a higher-order accuracy formula for the collocation method. As shown before, when the second-order accuracy formula is used, k is set as 3, the length of the knot sequence is N + 3, and MBρ,
Employing the backward Euler scheme for time integration, which is the first-order accurate in time,
(15) |
Based on FDM, the diffusion term is discretized by the second-order center difference as
The subscript i - 1/2 means value between grid points i - 1 and i, for example
(16) |
Inner boundary condition i = 1 is implemented by the second-order forward finite difference formula
(17) |
Outer boundary condition i = N,
(18) |
Finally, equations (16)–(18) constitute a matrix equation,
(19) |
in which MFDM is a tridiagonal matrix with one additional nonzero element (the third one in the first row), and
The matrix equation (14) derived by the collocation method and equation (19) by FDM can both be written like this:
(20) |
In general, M,
The predictor–corrector method generally suffers from a large numerical oscillation when coupled to theory-based transport models, such as the stiff transport model. The under-relaxation method and regula-falsi method can overcome this difficulty [24]. In this work, the under-relaxation method is applied if the predictor–corrector iteration can hardly converge. The diffusion term is calculated using χn+1, l like this:
(21) |
in which ϵ is the under-relaxation factor and superscript l stands for the l-th iteration at each time step.
Based on the FDM and collocation method mentioned in the previous section, two solvers (FDMS and BITS) are developed. Their results for different cases are compared in this section. Unless explicitly specified, BITS takes the same second-order accuracy (k = 3) in space as FDMS. For the first-order accuracy discretization in time and the second-order accuracy discretization in space, the discretization errors are linear over the time step and quadratic over the space step. Since iteration is usually needed to solve the nonlinear equation (20) for each time step, this brings additional errors. For simplicity, uniform ∆ρ in the minor radius is used from now on for both solvers; however, a non-uniform interval is also applicable.
For the sake of simplicity, density n and conduction coefficient χ are set as constant. In the steady state, the energy transport equation (1) is simplified as:
(22) |
This equation has the analytical solution
The numerical results of the temperature T by the FDM and collocation method are compared. The numerical error of T is calculated by the formula
(23) |
where Ti, analytic denotes the analytic solution of T at the i-th grid point.
As the grid points in ρ ∈ [0, 1] increase from 9 to 17, 33, 65, 129, the results of these two solvers both converge. The calculated errors by these two solvers are plotted in figure 1. As a cross-check, the results by Park (figure 2 of [24]) are also plotted, in which the cited second-order FDM result is marked as Park_FDM, and the cited fourth-order Interpolated Differential Operator (IDO) scheme result is marked as Park_IDO. Although the IDO method has a fourth-order accuracy, the highest power of the Hermite interpolation function used in equation (2) of [24] is 5, so we compare it to the BITS result with k = 6 (the highest power of polynomial is 5). Firstly, the results of FDMS developed in this work are as good as the results of the same method developed by Park. Secondly, the results of the fifth-order BITS (k = 6, marked as BITSk6) are as good as the results of the IDO method. Finally, the error of the second-order BITS (k = 3, marked as BITSk3) is about 1/4 to 1/5 of the error of the second-order FDMS.
Since density n, transport coefficient χ and source term S are independent of time and T in this case, then M and
Equation (1) is numerically solved by FDMS and BITS for two different χ models.
Initially, constant density and thermal diffusivity are used, n(ρ) = 1, χ(ρ) = 1, while the source takes the same form as in the previous case
The evolution of temperature T(ρ, t) is calculated by FDMS and BITS. The temperature profiles at the final time step by FDMS and BITS both converge as the number of grid points increases. In order to quantify the error, the N = 129 results are taken as standard, and the errors of other cases are calculated. For example, ρi in the N = 9 case is the same point ρ(i-1)×16+1 in the N = 129 case, so the calculated temperature in the N = 9 case Ti, 9 takes the same point result as in the N = 129 case T(i-1)×16+1, 129 as standard. Then the errors of N = 9 by FDMS and BITS are calculated separately by
(24) |
(25) |
The errors of other cases are calculated in a similar way. They are plotted at figure 2(a). The error of temperature by BITS is about 1/4–1/5 of that by FDMS. In order to reduce the influence of time discretization, this case is recalculated using a much smaller time step(dt decreases from 10-2 to 10-4), the result is plotted at figure 2(b). It is shown that the error by BITS is still about 1/4–1/5 of the error by FDMS for grid points N = 33, 65.
A stiff transport model is used to compare the performances of FDMS and BITS.
(26) |
Here κ represents the strength of the stiffness. χ0 = 1, κ = 10, γ = 1 and
The evolution of temperature T(ρ, t) is calculated by FDMS and BITS. The temperature profiles at the final time step by FDMS and BITS both converge as the number of grid points increases. The N = 129 results are taken as standard, and the errors of other cases are calculated as in the previous model. They are plotted in figure 3(a). The profile of χ at the last time step by FDMS and BITS is plotted in figure 3(b). The error by BITS is almost the same as by FDMS. In order to decrease the error from time discretization, the case is recalculated using a much smaller time step (dt decreases from 10-2 to 10-4). As shown in figure 3(c), the error of temperature by BITS decreases to about 1/4 to 1/5 of the error by FDMS for grid points N = 33, 65.
Collisional transport is common in tokamak plasmas. The coefficients of conductivity have been computed for completely ionized gases by [28]. The collisional transport coefficient varies as T-1/2. Here the transport coefficient is simplified as below:
(27) |
The total heating source is S = Sh + Qe, in which
The evolution of temperature T(ρ, t) is calculated by FDMS and BITS. The temperature profiles at the final time step by FDMS and BITS both converge as the number of grid points increases. The N = 129 results are taken as standard, and the errors of other cases are calculated as in the previous model. They are plotted in figure 4(a). The error by FDMS is a little better than by BITS at about 0.55–0.65 of BITS. In order to decrease the error from time discretization, the case is also recalculated using a much smaller time step (dt decreases from 10-3 to 10-5). As shown in figure 4(b), the error of temperature by BITS is almost the same as the error by FDMS.
In this subsection, the one-dimensional diffusion equations of density n, electron and ion temperature Te and Ti, and the poloidal magnetic field Bp are solved by these two methods. For simplicity, a low-β circular cross section plasma is considered, which means the geometry factors that appear at the surface average can be set as one. The equations of n, Te, Ti, and Bp are listed below:
(28) |
(29) |
(30) |
(31) |
in which
FDMS and BITS are also developed to solve these four equations. A similar algorithm is used in FDMS just like the previous subsection: a three-point, centered difference is used for the diffusion term, a two-point, upwind difference is used for the convection term, and a two-point difference is used for the time derivative. In BITS, the unknowns are arranged in this way:
(32) |
in which MDBρ is obtained by expanding all the elements of MBρ in equation (3) to a 4 × 4 diagonal matrix,
(33) |
(34) |
The derivative of variables is handled in a similar way. Then these four transport equations can be discretized into sparse linear systems and solved efficiently, just as we have done for the single transport equation case in the previous subsection.
The initial density, electron and ion temperature, and flux averaged toroidal current are set as below:
(35) |
(36) |
(37) |
\left\langle j_\phi R_0 / R\right\rangle\left(\rho, t_0\right)=\left(1-\rho^2\right)^2 | (38) |
The initial poloidal magnetic field Bp can be integrated from the flux averaged toroidal current as below:
\rho B_{\mathrm{p}}\left(\rho, t_0\right)=\mu_0 \int_0^\rho\left\langle j_\phi R_0 / R\right\rangle\left(\rho, t_0\right) \rho \mathrm{d} \rho | (39) |
The units of density, temperature and current are 1020 m-3, keV, and 106 A m-2, respectively. The sources are set as Sn = 0, Qe = Sh - Q∆, Qi = Q∆, in which Sh is auxiliary heat,
Then these four equation are solved by FDMS and BITS, using the predictor–corrector iteration with dt = 0.001 s, 100 steps. The magnetic equilibrium is fixed for simplicity. The error tolerance in iteration control is set roughly as 10-2, and only one corrector sub-step is needed for both FDMS and BITS.
The N = 129 results are taken as standard, and the errors of density n, electron temperature Te, ion temperature Ti, multiplication of minor radius and poloidal magnetic field ρBp are calculated as for the previous model. They are plotted in figure 5. The error of density by FDMS is about 0.6 of BITS for the N = 9 case. However, the performance of BITS is improved faster as the grid points increase. The error of density by BITS is about 1/2 that of FDMS for N = 65 case. The errors of Te and Ti by BITS are slightly smaller than the FDMS results for the N = 9 case. The performance of BITS is improved much faster than that of FDMS as the grid points increase. The error of temperature by BITS is about 1/5 that of FDMS for the N = 33 case and about 1/10 that of FDMS for the N = 65 case. This means the improvement of results in FDMS by increasing the grid points is deteriorated due to the coupling of transport equations. The density equation is weakly coupled to the Te equation, and the error of density by FDMS is slightly deteriorated. Te and Ti have much strong coupling with each other, and the errors of Te and Ti by FDMS are heavily deteriorated. The improvement of results in BITS by increasing grid points is still as good as in the single transport equation case (the slope of Te, Ti by BITS is as steep as the slope of density by FDMS). The diffusion coefficient of ρBp is set as independent of others, and much smaller than the density and thermal diffusion coefficients, so ρBp varies much more slowly. The error of ρBp by FDMS and BITS is very close for all the different grid points.
BITS is developed to solve 1-D diffusion equations for tokamak plasma based on the B-spline collocation method. It takes the same level of computation cost compared to the FDM which is widely used in existing tokamak transport codes.
The numerical results of a few cases are compared with the same order accuracy FDMS. The error of temperature by BITS is about 1/4 to 1/5 that of FDMS for the steady-state temperature transport equation case, in which no time discretization error or iteration error appears. For the single-variable transport equation, three different transport models are used. For the simple transport model (constant χ), the error by BITS is also about 1/4 to 1/5 that of FDMS. For the stiff and collision models, the performance of BITS becomes the same or even a little worse than that of FDMS. By using a smaller time step, the error by BITS decreases to 1/4 to 1/5 that of FDMS in the stiff case or about the same in the collision case. For coupled multi-variable transport equations (collisional transport coefficients are used), the error of variables by BITS is about the same as that of FDMS, if the variable (such as density or poloidal magnetic field) is not closely coupled with other unknown variables, while the error of variables by BITS is much smaller than FDMS, if the variables (such as electron and ion temperature) are strongly coupled with each other. The improvement in FDMS results by increasing grid points is deteriorated if they are strongly coupled with each other. The improvements in BITS results by increasing grid points are still as good as in the single transport equation case. Taking the coupled 1-D tokamak plasmas transport equations as an example, the maximum errors of n, Te, Ti, and ρBp by BITS using N = 33 are still about half of the maximum ones by FDMS using N = 65. This means BITS can reduce the computation cost by half (the computation cost of band diagonal linear equations is proportional to the number of grid points N), double the calculation speed and get more accurate results than FDMS. Since large computational resources are demanded when the theory-based turbulent transport model for tokamak plasmas is coupled to the 1-D transport solver to predict the time evolution of tokamak plasmas or its steady-state solution, the transport solver presented here can speed up the numerical analysis greatly.
The B-spline basis function can be stably calculated by the following recurrence relation.
B_{i, 1}= \begin{cases}1 & t_i \leqslant x<t_{i+1} \\ 0 & \text { otherwise }\end{cases} | (40) |
B_{i, k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i} B_{i, k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}} B_{i+1, k-1} | (41) |
Its derivative
(42) |
in which
(43) |
The first and second derivatives are listed below for k = 3, tl ≤ x < tl+1:
(44) |
(45) |
(46) |
(47) |
(48) |
(49) |
In order to give a vivid picture of the basis function B(x), let us specify a knot sequence as t1 = 1, ···, ti = i, ···, tn = n. Then for a given x which falls into the interval ti ≤ x < ti+1, the relevant k ≤ 3 order basis function will take the form as below.
(50) |
(51) |
(52) |
(53) |
(54) |
(55) |
Bi, 1, Bi-1, 2, Bi, 2, Bi-2, 3, Bi-1, 3, Bi, 3 have nonzero values in the interval [ti, ti+1), which are plotted in figure A1.
This work was supported by the National MCF Energy R&D Program of China (No. 2019YFE03040004), the Comprehensive Research Facility for Fusion Technology Program of China (No. 2018-000052-73-01-001228), and the National MCF Energy R&D Program of China (No. 2019YFE03060000).
[1] |
Shimizu K, Takizuka T and Sakasai A 1997 J. Nucl. Mater. 241–3 167 doi: 10.1016/S0022-3115(96)00503-X
|
[2] |
Connor J W et al 2004 Nucl. Fusion 44 R1 doi: 10.1088/0029-5515/44/4/R01
|
[3] |
Ida K and Rice J E 2014 Nucl. Fusion 54 045001 doi: 10.1088/0029-5515/54/4/045001
|
[4] |
Rosenbluth M N, Hazeltine R D and Hinton F L 1972 Phys. Fluids 15 116 doi: 10.1063/1.1693728
|
[5] |
Hinton F L and Hazeltine R D 1976 Rev. Mod. Phys. 48 239 doi: 10.1103/RevModPhys.48.239
|
[6] |
Pfeiffer W W et al 1980 ONETWO: A Computer Code for Modeling Plasma Transport in Tokamaks GA-A-16178 General Atomics, San Diego, CA
|
[7] |
Hawryluk R J 1981 An empirical approach to tokamak transport Physics of Plasmas Close to Thermonuclear Conditions (New York: Pergamon)
|
[8] |
Ongena J P H E, Evrard M and McCune D 1998 Fusion Technol. 33 181 doi: 10.13182/FST98-A11947009
|
[9] |
Pereverzev G V and Yushmanov P N 2002 ASTRA. Automated System for Transport Analysis in a Tokamak (IPP 5/98) (Garching: Max-Planck-Institut für Plasmaphysik)
|
[10] |
Coster D P et al 2010 IEEE Trans. Plasma Sci. 38 2085 doi: 10.1109/TPS.2010.2056707
|
[11] |
Gao X et al 2020 Nucl. Fusion 60 102001 doi: 10.1088/1741-4326/abaa91
|
[12] |
Zhu Y R et al 2022 Nucl. Fusion 62 076011 doi: 10.1088/1741-4326/ac65fb
|
[13] |
Zou Y P et al 2021 Chin. Phys. Lett. 38 045203 doi: 10.1088/0256-307X/38/4/045203
|
[14] |
Artaud J F et al 2010 Nucl. Fusion 50 043001 doi: 10.1088/0029-5515/50/4/043001
|
[15] |
Romanelli M et al 2014 Plasma Fusion Res. 9 3403023 doi: 10.1585/pfr.9.3403023
|
[16] |
Grierson B A et al 2018 Fusion Sci. Technol. 74 101 doi: 10.1080/15361055.2017.1398585
|
[17] |
Romanelli M et al 2020 Fusion Sci. Technol. 76 894 doi: 10.1080/15361055.2020.1819751
|
[18] |
Lee C Y et al 2021 Nucl. Fusion 61 096020 doi: 10.1088/1741-4326/ac1690
|
[19] |
Šesnić S et al 2018 IEEE Trans. Plasma Sci. 46 1027 doi: 10.1109/TPS.2018.2811858
|
[20] |
Soler J A et al 2020 J. Comput. Phys. 405 109093 doi: 10.1016/j.jcp.2019.109093
|
[21] |
Waltz R E et al 1997 Phys. Plasmas 4 2482 doi: 10.1063/1.872228
|
[22] |
Bateman G et al 1998 Phys. Plasmas 5 1793 doi: 10.1063/1.872848
|
[23] |
Staebler G M, Kinsey J E and Waltz R E 2007 Phys. Plasmas 14 055909 doi: 10.1063/1.2436852
|
[24] |
Park J M et al 2017 Comput. Phys. Commun. 214 1 doi: 10.1016/j.cpc.2016.12.018
|
[25] |
de Boor C 1978 A Practical Guide to Splines(New York: Springer)
|
[26] |
Dehghan M and Fakhar-Izadi F 2011 Math. Comput. Model. 53 1865 doi: 10.1016/j.mcm.2011.01.011
|
[27] |
Press W H et al 1986 Numerical Recipes: The Art of Scientific Computing 2nd edn (Cambridge: Cambridge University Press)
|
[28] |
Spitzer L J and Härm R 1953 Phys. Rev. 89 977 doi: 10.1103/PhysRev.89.977
|
1. | Zhu, Q., Zhang, S., Han, Q. et al. Theoretical kinetics with a one-dimensional fluid model and experimental investigation of coaxial XeCl excilamps. Journal of Physics D: Applied Physics, 2025, 58(7): 075201. DOI:10.1088/1361-6463/ad933d | |
2. | Qin, S., Wang, M., Du, J. et al. Effect of power supply parameters on discharge characteristics and sterilization efficiency of magnetically driven rotating gliding arc. Plasma Science and Technology, 2024, 26(9): 094006. DOI:10.1088/2058-6272/ad547d | |
3. |
Zhao, L.-Y., Zhang, J.-L., Zhang, K. et al. Experimental study on the inactivation of bioaerosols by dielectric barrier discharge combined with Ag-Cu/TiO2-CS composite photocatalyst | [介质阻挡放电协同 Ag-Cu/TiO2-CS 复合光催化剂灭活生物气溶胶的实验研究]. Zhongguo Huanjing Kexue/China Environmental Science, 2024, 44(5): 2777-2785.
![]() | |
4. | Chen, X., Li, M., Wang, W. et al. On the evolution and formation of discharge morphology in pulsed dielectric barrier discharge. Plasma Science and Technology, 2024, 26(4): 045403. DOI:10.1088/2058-6272/ad13e4 | |
5. | Li, Y., Lu, J., Gao, M. et al. Improvement of Vanadium Redox Flow Battery Efficiency Through Carbon Felt Electrodes Modification by Atmospheric Dielectric Barrier Discharge. Plasma Chemistry and Plasma Processing, 2023, 43(6): 1673-1684. DOI:10.1007/s11090-023-10365-4 |