BITS: an efficient transport solver based on a collocation method with B-spline basis

Xiaotao XIAO; Shaojie WANG; Lei YE; Zongliang DAI; Chengkang PAN; Qilong REN

doi:10.1088/2058-6272/ac909f

Plasma Science and Technology > 2023 > 25(2): 025103. > DOI: 10.1088/2058-6272/ac909f

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

Citation:

PDF (525 KB)

BITS: an efficient transport solver based on a collocation method with B-spline basis

1.
Institute of Plasma Physics, Chinese Academy of Science, Hefei 230031, People's Republic of China
2.
Department of Modern Physics, University of Science and Technology of China, Hefei 230026, People's Republic of China

More Information

Corresponding author:
Xiaotao XIAO, E-mail: xtxiao@ipp.ac.cn
Received Date: May 19, 2022
Revised Date: July 23, 2022
Accepted Date: September 07, 2022
Available Online: December 05, 2023
Published Date: January 05, 2023

Graphical Abstract

Abstract

Abstract

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.
- plasma simulation,
- transport equation,
- collocation,
- B-spline,
- tokamak

FullText(HTML)

1. Introduction

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.

2. Description of scheme

The 1-D energy diffusion equation for magnetically confined toroidal plasma in an axisymmetric tokamak can be simplified in cylindrical coordinates as

$\frac{3}{2} n \frac{\partial T}{\partial t} - \frac{1}{ρ} \frac{\partial}{\partial ρ} [ρ (n χ \frac{\partial T}{\partial ρ})] = 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 ${T^{'}|}_{ρ = 0} = 0,$ ${T|}_{ρ = 1} = T_{b}$ , and initial condition T(ρ, t = 0) = T₀(ρ), this equation is well posed. The radial coordinate ρ is uniformly divided into grid points ρ₁ = 0, ···, ρ_i, ···, ρ_N = 1 for simplicity, i.e. ∆ρ = ρ_i+1 - ρ_i = (ρ_N - ρ₁)/(N - 1), i = 1, ···, N - 1. A non-uniform interval is also applicable. T at time step n is marked as Tⁿ.

2.1 Collocation method

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 B_j(x) is equal to the given value y_i at points x_i, i.e. ∑_jα_jB_j(x_i) = y_i. When solving the differential equation, the constraint is that the combination of coefficients α_j, basis functions B_j(x) and their derivatives $B_{j}^{'}$ , $B_{j}^{'}$ and so on is equal to zero, i.e. $c_{0} \sum_{j} α_{j} B_{j} (x_{i}) + c_{1} \sum_{j} α_{j} B_{j}^{'} (x_{i}) + \cdot\cdot\cdot = 0$ . Here c₀, c₁ stand for the factors appearing in the differential equation $c_{0} y + c_{1} y^{'} + \cdot\cdot\cdot = 0$ . This is the basic idea of the interpolation method (collocation method) to solve differential equations. This method has been used to solve ordinary differential equations (chapter 15 of [25]), and Ostrovsky equations [26].

Based on collocation, T is approximated by a linear combination of basis function B_{i, k}(ρ). Here the second-order accuracy is used, i.e. T(ρ) = Σ_iα_iB_{i, 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 t₁, ···, t_N+k is generated, then B_{i, k}(ρ) is totally determined. In this work, the knot sequence is generated as follows:

$\{\begin{matrix} t_{1} = t_{2} = t_{3} = ρ_{1}, \\ t_{i} = (ρ_{i - 2} + ρ_{i - 1}) / 2, & i = 4, \cdot\cdot\cdot, N \\ t_{N + 1} = t_{N + 2} = t_{N + 3} = ρ_{N} . \end{matrix}$

(2)

Then temperature T at grid points $\vec{T} = (T_{1}, \cdot\cdot\cdot, T_{N})$ can be represented by coefficient $\vec{α}$ as below (it is also the matrix equation solved in a typical interpolation problem):

$\begin{matrix} (\begin{matrix} T_{1} \\ ⋮ \\ T_{N} \end{matrix}) & = & (\begin{matrix} B_{1} (ρ_{1}) & \dots & B_{N} (ρ_{1}) \\ ⋮ & ⋱ & ⋮ \\ B_{1} (ρ_{N}) & \dots & B_{N} (ρ_{N}) \end{matrix}) \cdot (\begin{matrix} α_{1} \\ ⋮ \\ α_{N} \end{matrix}) \\ = & M_{B ρ} \cdot \vec{α} \end{matrix}$

(3)

in which M_Bρ is a tridiagonal matrix, with only two additional nonzero elements m_{1, 3}, m_{N, N-2}. Here element M_Bρ(i, j) is written as m_{i, j} = B_j(ρ_i).

$M_{B ρ} = [\begin{matrix} m_{1, 1} & m_{1, 2} & m_{1, 3} \\ m_{2, 1} & m_{2, 2} & m_{2, 3} \\ ⋱ & ⋱ & ⋱ \\ m_{N, N - 2} & m_{N, N - 1} & m_{N, N} \end{matrix}]$

(4)

As temperature T varies with time, so does the coefficient. We have $T^{n} = M_{B ρ} \cdot {\vec{α}}^{n}$ at time step n, $T^{n + 1} = M_{B ρ} \cdot {\vec{α}}^{n + 1}$ at time step n + 1. Hence, the time-varying term discretized by the backward finite difference is represented as:

$\begin{matrix} \frac{3}{2} n \frac{\partial T^{n + 1}}{\partial t} = \frac{3}{2} n \frac{T^{n + 1} - T^{n}}{∆ t} \\ = \frac{3 n}{2 ∆ t} M_{B ρ} \cdot ({\vec{α}}^{n + 1} - {\vec{α}}^{n}) \end{matrix}$

(5)

${\vec{T}}^{'} = (T_{1}^{'}, \cdot\cdot\cdot, T_{N}^{'})$ can be represented as

${\vec{T}}^{'} = \frac{d M_{B ρ}}{d ρ} \cdot \vec{α} = M_{B ρ'} \cdot \vec{α},$

(6)

in which

$M_{B ρ'} = [\begin{matrix} m_{1, 1}^{'} & m_{1, 2}^{'} & m_{1, 3}^{'} \\ m_{2, 1}^{'} & m_{2, 2}^{'} & m_{2, 3}^{'} \\ ⋱ & ⋱ & ⋱ \\ m_{N, N - 2}^{'} & m_{N, N - 1}^{'} & m_{N, N}^{'} \end{matrix}]$

(7)

${\vec{T}}^{''} = (T_{1}^{''}, \cdot\cdot\cdot, T_{N}^{''})$ can be represented as

${\vec{T}}^{''} = \frac{d^{2} M_{B ρ}}{d ρ^{2}} \cdot \vec{α} = M_{B ρ''} \cdot \vec{α},$

(8)

in which

$M_{B ρ''} = [\begin{matrix} m_{1, 1}^{''} & m_{1, 2}^{''} & m_{1, 3}^{''} \\ m_{2, 1}^{''} & m_{2, 2}^{''} & m_{2, 3}^{''} \\ ⋱ & ⋱ & ⋱ \\ m_{N, N - 2}^{''} & m_{N, N - 1}^{''} & m_{N, N}^{''} \end{matrix}]$

(9)

Here M_Bρ, $M_{B ρ'}$ , M_Bρ'' are all tridiagonal matrices, and only two additional elements are nonzero: the third one in the first row and the third from last one in the last row. Since grid points ρ₁, ···, ρ_N and knot sequence t₁, ···, t_N+3 remain unchanged after initialization, all the elements of matrix M_Bρ, $M_{B ρ'}$ , M_Bρ'' can be calculated once and used repeatedly.

For interior points i = 2, ···, N - 1, the equation is discretized as

$\begin{matrix} (\frac{3 n}{2 ∆ t} M_{B ρ} - n χ^{n + 1} M_{B ρ''} - \frac{1}{ρ} \frac{\partial ρ n χ^{n + 1}}{\partial ρ} M_{B ρ'}) \cdot {\vec{α}}^{n + 1} \\ = \frac{3 n}{2 ∆ t} M_{B ρ} \cdot {\vec{α}}^{n} + S \end{matrix}$

(10)

For the inner boundary condition ${T^{'}|}_{ρ = 0} = 0$ ,

$α_{1} B_{1}^{'} (ρ_{1}) + α_{2} B_{2}^{'} (ρ_{1}) + α_{3} B_{3}^{'} (ρ_{1}) = 0$

(11)

For the outer boundary condition ${T|}_{ρ = 1} = T_{b}$ ,

$\begin{matrix} α_{N - 2} B_{N - 2} (ρ_{N}) + α_{N - 1} B_{N - 1} (ρ_{N}) \\ + α_{N} B_{N} (ρ_{N}) = T_{b} \end{matrix}$

(12)

It is easy to construct a general boundary condition in a similar way. For example, $(β_{1} T + β_{2} T^{'}) |_{ρ = 1} = T_{b}$ can be constructed as

$\begin{matrix} α_{N - 2} [β_{1} B_{N - 2} (ρ_{N}) + β_{2} B_{N - 2}^{'} (ρ_{N})] \\ + α_{N - 1} [β_{1} B_{N - 1} (ρ_{N}) + β_{2} B_{N - 1}^{'} (ρ_{N})] \\ + α_{N} [β_{1} B_{N} (ρ_{N}) + β_{2} B_{N}^{'} (ρ_{N})] = T_{b} \end{matrix}$

(13)

Finally, equations (10)–(12) constitute a matrix equation,

$M_{CM} \cdot {\vec{α}}^{n + 1} = {\vec{b}}_{CM}$

(14)

in which M_CM 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); ${\vec{α}}^{n + 1}$ and ${\vec{b}}_{CM}$ are column vectors. It can also be solved efficiently by lower–upper decomposition [27].

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 M_Bρ, $M_{B ρ'}$ , M_Bρ'' have three nonzero elements in each row. If we want to use the fourth-order accuracy formula, we simply set k as 5, generate a proper knot sequence with a length of N + 5, and matrices M_Bρ, $M_{B ρ'}$ , M_Bρ'' have five nonzero elements in each row, while all other things remain the same.

2.2 FDM

Employing the backward Euler scheme for time integration, which is the first-order accurate in time,

$\frac{3}{2} n {\frac{\partial T^{n + 1}}{\partial t}|}_{i} = \frac{3}{2} n \frac{T_{i}^{n + 1} - T_{i}^{n}}{∆ t}$

(15)

Based on FDM, the diffusion term is discretized by the second-order center difference as

$\frac{1}{ρ_{i}} \frac{1}{∆ ρ} [{(ρ n χ)}_{i + 1 / 2} \frac{T_{i + 1} - T_{i}}{∆ ρ} - {(ρ n χ)}_{i - 1 / 2} \frac{T_{i} - T_{i - 1}}{∆ ρ}]$

The subscript i - 1/2 means value between grid points i - 1 and i, for example ${(ρ n χ)}_{1 + 1 / 2} = [{(ρ n χ)}_{i} + {(ρ n χ)}_{i + 1}] / 2$ . The resulting equation for interior points i = 2, ···, N - 1 is written as follows:

$\begin{matrix} - \frac{1}{ρ_{i}} \frac{1}{∆ ρ^{2}} {(ρ n χ)}_{i - 1 / 2} T_{i - 1}^{n + 1} \\ + \{\frac{3 n}{2 ∆ t} + \frac{1}{ρ_{i}} \frac{1}{∆ ρ^{2}} [{(ρ n χ)}_{i + 1 / 2} - {(ρ n χ)}_{i - 1 / 2}]\} T_{i}^{n + 1} \\ - \frac{1}{ρ_{i}} \frac{1}{∆ ρ^{2}} {(ρ n χ)}_{i + 1 / 2} T_{i + 1}^{n + 1} \\ = \frac{3}{2} n \frac{T_{i}^{n}}{∆ t} + S_{i} \end{matrix}$

(16)

Inner boundary condition i = 1 is implemented by the second-order forward finite difference formula

$T^{'} (ρ_{1}) = \frac{- 1.5 T_{1} + 2 T_{2} - 0.5 T_{3}}{∆ ρ} = 0$

(17)

Outer boundary condition i = N,

$T_{N} = T (ρ_{N}, t)$

(18)

Finally, equations (16)–(18) constitute a matrix equation,

$M_{FDM} \cdot {\vec{T}}^{n + 1} = {\vec{b}}_{FDM}$

(19)

in which M_FDM is a tridiagonal matrix with one additional nonzero element (the third one in the first row), and ${\vec{T}}^{n + 1}$ and ${\vec{b}}_{FDM}$ are column vectors. It can be solved efficiently by lower–upper decomposition [27].

2.3 Iteration method

The matrix equation (14) derived by the collocation method and equation (19) by FDM can both be written like this:

$M (T^{n + 1}) \cdot {\vec{T}}^{n + 1} = \vec{b} (T^{n + 1})$

(20)

In general, M, $\vec{b}$ depend upon Tⁿ⁺¹ (as well as $T^{' n + 1}$ ). Hence, it is nonlinear with regard to Tⁿ⁺¹. This nonlinearity is usually handled by a predictor–corrector method. In the predictor sub-step, M, $\vec{b}$ are evaluated by M(Tⁿ), $\vec{b} (T^{n})$ , then T^{n+1, 1} is solved by equation (20). In the corrector sub-step, M, $\vec{b}$ are evaluated by M(T^{n+1, 1}), $\vec{b} (T^{n + 1, 1})$ using the predicted T^{n+1, 1}, then T^{n+1, 2} is obtained by solving equation (20). If the predicted value T^{n+1, 1} is close enough to the corrected value T^{n+1, 2}, it then proceeds to the next time step. Otherwise the corrected value is used as a new predicted value and the corrector sub-step is iterated until convergence is achieved or until the maximum number of iterations is exceeded.

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:

$χ^{n + 1, l} = \{\begin{matrix} χ (T^{' n}) & l = 1 \\ ϵ χ (T^{' n + 1, l}) + (1 - ϵ) χ^{n + 1, l - 1} & l > 1 \end{matrix}$

(21)

in which ϵ is the under-relaxation factor and superscript l stands for the l-th iteration at each time step.

3. Results

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.

3.1 Steady-state transport equation

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:

$n χ T^{'} + n χ T^{'} / ρ + S = 0$

(22)

This equation has the analytical solution $T = \frac{S_{0}}{n χ} (\exp (1 - ρ^{2}) - 1)$ for $S = 4 S_{0} (1 - ρ^{2}) \exp (1 - ρ^{2})$ , while the boundary conditions are applied as ${T^{'}|}_{ρ = 0} = 0$ and ${T|}_{ρ = 1} = 0$ . S₀ = 1 is used in the work.

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

$σ_{T} = \frac{Σ_{i = 1}^{N} |T_{i} - T_{i, analytic}|}{Σ_{i = 1}^{N} |T_{i, analytic}|}$

(23)

where T_{i, 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.

Figure 1. Error of T by different solvers. For the second-order scheme, FDMS is as good as Park_FDM, and the error of BITS(BITSk3) is about 1/4 to 1/5 that of FDMS. The fifth-order BITS with k = 6 (BITSk6) is as good as Park_IDO.

DownLoad: Full-Size Img PowerPoint

Since density n, transport coefficient χ and source term S are independent of time and T in this case, then M and $\vec{b}$ in equation (20) are both independent of Tⁿ⁺¹. This equation is simplified to a linear one, and the error introduced by iteration disappears. There is no time-varying term, so the error introduced by time discretization also disappears. Only the error from space discretization is retained. Hence, this case shows the second-order BITS decrease the error from space discretization to 1/4 to 1/5 that of the same order FDMS.

3.2 Single-variable transport equation

Equation (1) is numerically solved by FDMS and BITS for two different χ models.

3.2.1 χ = 1 model

Initially, constant density and thermal diffusivity are used, n(ρ) = 1, χ(ρ) = 1, while the source takes the same form as in the previous case $S = 4 S_{0} (1 - ρ^{2}) \exp (1 - ρ^{2})$ . The initial condition is T(ρ, t = 0) =0. T approaches the steady state $T = \frac{S_{0}}{n χ} (\exp (1 - ρ^{2}) - 1)$ . The time step dt = 0.01, total steps 100, and different grid points 9, 17, 33, 65, 129 are used. The error tolerance in iteration control is set as 10^-4, and only one corrector sub-step is needed for both FDMS and BITS.

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 T_{i, 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

$σ_{9, FDMS} = \frac{Σ_{i = 1}^{9} |T_{i, 9, FDMS} - T_{(i - 1) \times 16 + 1, 129, FDMS}|}{Σ_{i = 1}^{9} |T_{(i - 1) \times 16 + 1, 129, FDMS}|}$

(24)

$σ_{9, BITS} = \frac{Σ_{i = 1}^{9} |T_{i, 9, BITS} - T_{(i - 1) \times 16 + 1, 129, BITS}|}{Σ_{i = 1}^{9} |T_{(i - 1) \times 16 + 1, 129, BITS}|}$

(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.

Figure 2. The error by BITS is about 1/4 to 1/5 of the error by FDMS for dt = 10^-2 as shown in (a). For dt = 10^-4, the error by BITS is still about 1/4 to 1/5 of the error by FDMS for grid points N = 33, 65 as shown in (b).

DownLoad: Full-Size Img PowerPoint

3.2.2 Stiff model

A stiff transport model is used to compare the performances of FDMS and BITS.

$χ = \{\begin{matrix} χ_{0} + κ {(|T^{'}| - |T_{c}^{'}|)}^{γ}, & |T^{'}| > |T_{c}^{'}| \\ χ_{0}, & |T^{'}| < |T_{c}^{'}| \end{matrix}$

(26)

Here κ represents the strength of the stiffness. χ₀ = 1, κ = 10, γ = 1 and $T_{c}^{'} = 0.5$ are used. The initial condition is T = 0 and $S = 4 (1 - ρ^{2}) \exp (1 - ρ^{2})$ . The density is still set as constant n(ρ) = 1. FDMS and BITS using the predictor–corrector iteration both encounter convergence difficulty. The under-relaxation iteration is applied to overcome this difficulty. The under-relaxation factor ϵ is set as 0.3, the error tolerance in iteration control is set as 10^-4, and the iteration times in FDMS and BITS are similar. The same time step dt = 0.01, total steps 100, and different grid points 9, 17, 33, 65, 129 are used.

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.

Figure 3. The error by BITS is the same as FDMS for dt = 10^-2 as shown in (a). The χ at the last time step is shown in (b). For dt = 10^-4, the error by BITS is about 1/4 to 1/5 of the error by FDMS for grid points N = 33, 65 as shown in (c).

DownLoad: Full-Size Img PowerPoint

3.2.3 Collision model

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:

$χ = T^{- 1 / 2}$

(27)

The total heating source is S = S_h + Q_e, in which $S_{h} = 0.1 \exp (- \frac{1}{2} \frac{ρ^{2}}{{0.2}^{2}})$ is auxiliary heating and $Q_{e} = \frac{T - 0.1}{T^{3 / 2}}$ is the electron–ion energy exchange term. The initial condition is set as $T (ρ, t_{0}) = (3 - 0.1) \times {(1 - ρ^{2})}^{2} + 0.1$ . The time step dt = 0.001, total steps 100, and different grid points 9, 17, 33, 65, 129 are used. The error tolerance in iteration control is set as 10^-3, and the predictor–corrector iteration is used in both FDMS and BITS.

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.

Figure 4. The error by FDMS is about 0.55–0.65 of BITS for dt = 10^-3 as shown in (a). For dt = 10^-5, the error by BITS is almost the same as the error by FDMS as shown in (b).

DownLoad: Full-Size Img PowerPoint

3.3 Multi-variable transport equations

In this subsection, the one-dimensional diffusion equations of density n, electron and ion temperature T_e and T_i, and the poloidal magnetic field B_p 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, T_e, T_i, and B_p are listed below:

$\frac{\partial n}{\partial t} + \frac{\partial}{\partial ρ} (ρ Γ_{i}) = S_{n}$

(28)

$\frac{3}{2} T_{e} \frac{\partial n}{\partial t} + \frac{3}{2} n \frac{\partial T_{e}}{\partial t} + \frac{1}{ρ} \frac{\partial}{\partial ρ} [ρ (q_{e} + \frac{5}{2} T_{e} Γ_{e})] = Q_{e}$

(29)

$\frac{3}{2} T_{i} \frac{\partial n}{\partial t} + \frac{3}{2} n \frac{\partial T_{i}}{\partial t} + \frac{1}{ρ} \frac{\partial}{\partial ρ} [ρ (q_{i} + \frac{5}{2} T_{i} Γ_{i})] = Q_{i}$

(30)

$\frac{1}{ρ^{2}} \frac{1}{\partial t} (ρ B_{p}) - \frac{1}{ρ} \frac{\partial}{\partial ρ} [\frac{η}{μ_{0} ρ} \frac{\partial}{\partial ρ} (ρ B_{p})] = 0$

(31)

in which $Γ_{e} = - D \partial_{ρ} n_{e}$ , $Γ_{i} = - D \partial_{ρ} n_{i}$ are particle fluxes of electrons and ions; $q_{e} = - n χ_{e} \partial_{ρ} T_{e}$ , $q_{i} = - n χ_{i} \partial_{ρ} T_{i}$ are heat fluxes of electrons and ions; S_n, Q_e, Q_i are the sources of density, electron and ion thermal energy, respectively; and η is resistivity.

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: $\vec{u} (ρ_{1}), \cdot\cdot\cdot, \vec{u} (ρ_{N})$ , in which $\vec{u} (ρ_{i}) = (n (ρ_{i}), T_{e} (ρ_{i}), T_{i} (ρ_{i}), ρ B_{p} (ρ_{i}))$ , and then the unknowns can be represented as a product of sparse matrix M_DBρ and coefficient ${\vec{α}}_{V}$ :

$\begin{matrix} (\begin{matrix} \vec{u} (ρ_{1}) \\ ⋮ \\ \vec{u} (ρ_{N}) \end{matrix}) & = & (\begin{matrix} D_{B 1} (ρ_{1}) & \dots & D_{BN} (ρ_{1}) \\ ⋮ & ⋱ & ⋮ \\ D_{B 1} (ρ_{N}) & \dots & D_{BN} (ρ_{N}) \end{matrix}) \cdot (\begin{matrix} {\vec{α}}_{1} \\ ⋮ \\ {\vec{α}}_{N} \end{matrix}) \\ = & M_{DB ρ} \cdot {\vec{α}}_{V} \end{matrix}$

(32)

in which M_DBρ is obtained by expanding all the elements of M_Bρ in equation (3) to a 4 × 4 diagonal matrix,

$D_{Bj} (ρ_{i}) = [\begin{matrix} B_{j} (ρ_{i}) \\ B_{j} (ρ_{i}) \\ B_{j} (ρ_{i}) \\ B_{j} (ρ_{i}) \end{matrix}]$

(33)

${\vec{α}}_{V}$ is obtained by expanding all the elements of $\vec{α}$ in equation (3) to a vector,

${\vec{α}}_{V} = (\begin{matrix} {\vec{α}}_{1} \\ ⋮ \\ {\vec{α}}_{N} \end{matrix}), {\vec{α}}_{i} = (\begin{matrix} α_{i, n} \\ α_{i, T_{e}} \\ α_{i, T_{i}} \\ α_{i, ρ B_{p}} \end{matrix})$

(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:

$n (ρ, t_{0}) = (1 - 0.1) \times {(1 - ρ^{2})}^{2} + 0.1$

(35)

$T_{e} (ρ, t_{0}) = (3 - 0.1) \times {(1 - ρ^{2})}^{2} + 0.1$

(36)

$T_{i} (ρ, t_{0}) = (3 - 0.1) \times {(1 - ρ^{2})}^{2} + 0.1$

(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 B_p 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 10²⁰ m^-3, keV, and 10⁶ A m^-2, respectively. The sources are set as S_n = 0, Q_e = S_h - Q_∆, Q_i = Q_∆, in which S_h is auxiliary heat, $Q_{∆} = c_{∆} (T_{e} - T_{i}) = 3 m_{e} n (T_{e} - T_{i}) / (m_{i} τ_{e})$ is the electron–ion energy exchange term, and τ_e is the electron collision time. For singly charged ions, $τ_{e} = 6.4 \times 10^{14} T_{e}^{3 / 2} / n s$ , with T_e in keV. The transport coefficients take a simplified collisional form as $χ_{e} = 1.0 / T_{e}^{- 1 / 2} m^{2} s^{- 1}$ , $χ_{i} = 1.0 / T_{i}^{- 1 / 2} m^{2} s^{- 1}$ , $D = 0.25 / T_{e}^{- 1 / 2}$ m² s^-1, T_e and T_i in keV. η/μ₀ = 0.01 m s^-1. S_h is set as $S_{h} (ρ) = 0.1 \times \exp (- ρ^{2} / (2 \times {0.2}^{2})) M W$ .

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 T_e, ion temperature T_i, multiplication of minor radius and poloidal magnetic field ρB_p 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 T_e and T_i 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 T_e equation, and the error of density by FDMS is slightly deteriorated. T_e and T_i have much strong coupling with each other, and the errors of T_e and T_i 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 T_e, T_i by BITS is as steep as the slope of density by FDMS). The diffusion coefficient of ρB_p is set as independent of others, and much smaller than the density and thermal diffusion coefficients, so ρB_p varies much more slowly. The error of ρB_p by FDMS and BITS is very close for all the different grid points.

Figure 5. Error of n, T_e, T_i, ρB_p at the final time step by FDMS and BITS.

DownLoad: Full-Size Img PowerPoint

4. Summary

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, T_e, T_i, and ρB_p 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.

Appendix

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

$D^{j} (\sum_{i} α_{i} B_{i, k}) = \sum_{i} α_{i}^{(j + 1)} B_{i, k - j}$

(42)

in which

$α_{r}^{(j + 1)} = \{\begin{matrix} α_{r}, & j = 0 \\ \frac{α_{r}^{(j)} - α_{r - 1}^{(j)}}{(t_{r + k - j} - t_{r}) / (k - j)}, & j > 0 \end{matrix}$

(43)

The first and second derivatives are listed below for k = 3, t_l ≤ x < t_l+1:

$B_{l - 2, 3}^{'} (x) = - 2 \frac{B_{l - 1, 2} (x)}{(t_{l + 1} - t_{l - 1})}$

(44)

$B_{l - 1, 3}^{'} (x) = 2 \frac{B_{l - 1, 2} (x)}{(t_{l + 1} - t_{l - 1})} - 2 \frac{B_{l, 2} (x)}{(t_{l + 2} - t_{l})}$

(45)

$B_{l, 3}^{'} (x) = 2 \frac{B_{l, 2} (x)}{(t_{l + 2} - t_{l})}$

(46)

$B_{l - 2, 3}^{'} (x) = - 2 \frac{B_{l - 1, 2}^{'} (x)}{(t_{l + 1} - t_{l - 1})}$

(47)

$B_{l - 1, 3}^{'} (x) = 2 \frac{B_{l - 1, 2}^{'} (x)}{(t_{l + 1} - t_{l - 1})} - 2 \frac{B_{l, 2}^{'} (x)}{(t_{l + 2} - t_{l})}$

(48)

$B_{l, 3}^{'} (x) = 2 \frac{B_{l, 2}^{'} (x)}{(t_{l + 2} - t_{l})}$

(49)

In order to give a vivid picture of the basis function B(x), let us specify a knot sequence as t₁ = 1, ···, t_i = i, ···, t_n = n. Then for a given x which falls into the interval t_i ≤ x < t_i+1, the relevant k ≤ 3 order basis function will take the form as below.

$B_{i, 1} (x) = 1, i \leq x < i + 1$

(50)

$\begin{matrix} B_{i, 2} (x) & = & \frac{x - t_{i}}{t_{i + 1} - t_{i}} B_{i, 1} (x) + \frac{t_{i + 2} - x}{t_{i + 2} - t_{i + 1}} B_{i + 1, 1} (x) \\ = & x - i \end{matrix}$

(51)

$\begin{matrix} B_{i - 1, 2} (x) & = & \frac{x - t_{i - 1}}{t_{i} - t_{i - 1}} B_{i - 1, 1} (x) + \frac{t_{i + 1} - x}{t_{i + 1} - t_{i}} B_{i, 1} (x) \\ = & i + 1 - x \end{matrix}$

(52)

$\begin{matrix} B_{i, 3} (x) & = & \frac{x - t_{i}}{t_{i + 2} - t_{i}} B_{i, 2} (x) + \frac{t_{i + 3} - x}{t_{i + 3} - t_{i + 1}} B_{i + 1, 2} (x) \\ = & \frac{{(x - i)}^{2}}{2} \end{matrix}$

(53)

$\begin{matrix} B_{i - 1, 3} (x) & = & \frac{x - t_{i - 1}}{t_{i + 1} - t_{i - 1}} B_{i - 1, 2} (x) + \frac{t_{i + 2} - x}{t_{i + 2} - t_{i}} B_{i, 2} (x) \\ = & 1 - \frac{{(x - i)}^{2}}{2} - \frac{{(i + 1 - x)}^{2}}{2} \end{matrix}$

(54)

$\begin{matrix} B_{i - 2, 3} (x) & = & \frac{x - t_{i - 2}}{t_{i} - t_{i - 2}} B_{i - 2, 2} (x) + \frac{t_{i + 1} - x}{t_{i + 1} - t_{i - 1}} B_{i - 1, 2} (x) \\ = & \frac{{(i + 1 - x)}^{2}}{2} \end{matrix}$

(55)

B_{i, 1}, B_{i-1, 2}, B_{i, 2}, B_{i-2, 3}, B_{i-1, 3}, B_{i, 3} have nonzero values in the interval [t_i, t_i+1), which are plotted in figure A1.

A1. B(x).

DownLoad: Full-Size Img PowerPoint

Acknowledgments

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).

References (28)

References

[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

Cited By

Cited by

Periodical cited type(1)

1.	Ye, L., Luo, Z., Xiao, X. et al. Electron thermal internal transport barriers triggered by the effect of ion shielding. Nuclear Fusion, 2022, 62(12): 124002. DOI:10.1088/1741-4326/ac9c18