2.1. Problem definition

In this subsection, the assumptions of the computation are organised and the self-consistent computation is formulated in the form of an optimisation problem.

At a certain moment $t=t_0$, we assume that a consistent distribution of eddy currents and an MHD equilibrium have been obtained. When any temporal changes occur in the plasma, how can we obtain the correct solution at the next time step, $t=t_1$, using the VMEC code? For a clearer problem definition, an overview of the VMEC calculation is provided first.

The basic input variables for VMEC include the radial profiles of toroidal current and pressure, along with the total toroidal magnetic flux enclosed by the plasma boundary. Additionally, when determining the plasma shape through equilibrium calculations, the distribution of vacuum magnetic field (the magnetic field in the entire computational domain produced by all currents such as eddy, those flowing at coils, etc. except the plasma current) must be specified. Alternative options for input variables exist but are not used in this study. The calculation performed by VMEC determines the shape of the magnetic surfaces and the magnetic field distribution within the plasma so that the given radial profiles of toroidal current and pressure on the nested magnetic surfaces satisfy MHD equilibrium. The toroidal magnetic flux enclosed by the magnetic surfaces is used as the radial coordinate (magnetic surface label), and the input total toroidal magnetic flux corresponds to that value at the plasma boundary. The shape of the plasma boundary is determined such that force balance is achieved under the total magnetic field distribution, which is the sum of the input vacuum magnetic field and the magnetic field produced by the plasma current. In this study, the radial profile of pressure, the total toroidal magnetic flux and the components of the vacuum magnetic field produced by the coil currents are fixed at their values of $t_0$. Based on this setting, the problem solved in this paper can be rephrased as follows: we wish to determine the evolution of eddy currents consistent with the given temporal change of the radial profile of toroidal plasma current at $t_1$ under the condition that the plasma maintains MHD equilibrium. Note that, even though the radial profile of toroidal current is given, the plasma current distribution in the real space is not fixed, which is affected by the magnetic field produced by the eddy currents at $t_1$ and changes according to the deformation and displacements in the magnetic surfaces.

On the other hand, the KEDDY3D code, which computes the evolution of the eddy current distribution, requires the temporal variation of flux linkage on the conductor as an input. The temporal variation of the magnetic flux due to the plasma current is computed based on the equilibria output by VMEC at $t = t_0$ and $t_1$. The latter introduces an implicit relation between the VMEC and KEDDY3D codes, requiring the iterative approach.

Figure 1 illustrates the basic flow chart of the evolution calculation from $t = t_0$ to $t_1$. The magnetic field produced by the eddy currents, $\boldsymbol {B}^\text {eddy}_{t_1}$, is treated as an unknown variable. At the beginning of the calculation, an initial guess of $\boldsymbol {B}^\text {eddy}_{t_1}$ is provided as ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$ (for example, ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1} = \boldsymbol {B}^\text {eddy}_{t_0}$). Additionally, other input variables for equilibrium calculation at $t_1$ are also supplied (in the figure, only the coil magnetic field is explicitly shown). Using ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$, we obtain the plasma current distribution, $\boldsymbol {J}^\text {plasma}$, in real coordinates by VMEC calculation. Then, through numerical integration based on the Biot–Savart law, the magnetic flux linkage on the VV, $\varPsi ^\text {plasma}_{t_1}$, is calculated. The temporal difference of $\varPsi ^\text {plasma}$ allows for the temporal evolution of the eddy current distribution to be determined by KEDDY3D, finally yielding the magnetic field produced by the eddy currents, $\boldsymbol {B}^\text {eddy}_{t_1}$, as the output. If ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$ and $\boldsymbol {B}^\text {eddy}_{t_1}$ are sufficiently close, the calculation at t=t1 is regarded as being converged. However, if ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$ and $\boldsymbol {B}^\text {eddy}_{t_1}$ differ, the provisional magnetic field ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$ is updated and the calculation is repeated. The main focus of this paper is how to update ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$.

Figure 1. Flow of the evolution calculation to obtain the magnetic field produced by consistent eddy currents at $t=t_1$. Here, $\boldsymbol {J}^{\text {plasma}}$ and $\boldsymbol {J}^{\text {eddy}}$ represent the plasma current and the eddy currents, $\boldsymbol {B}^{\text {eddy}}$ denotes the magnetic field produced by $\boldsymbol {J}^{\text {eddy}}$, $\varPsi ^{\text {plasma}}$ denotes the magnetic flux produced by $\boldsymbol {J}^{\text {plasma}}$ and $\boldsymbol {B}^{\text {coil}}$ denotes the magnetic field produced by the coil current. The subscript of each variable indicates the value at the time step. The processes enclosed by the dashed line correspond to the calculations represented by function $\boldsymbol {f}$ in the text.

For updating ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$, the simplest approach is to replace the provisional magnetic field ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$ by the output $\boldsymbol {B}^\text {eddy}_{t_1}$ from the previous iteration. Before introducing the proposed iteration scheme, we consider the consequence of such a simple update scheme. This corresponds to providing $\boldsymbol {B}^\text {eddy}_{t_1}$ directly through the blue ‘Update’ box in figure 1 without any modification, allowing it to become the next ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1}$. When using this simple update method, the calculation often diverges instead of converging, due to a nature of the VMEC's calculation. As an example, consider applying a magnetic perturbation, referred to as a ‘kick’, directed downwards to a vertically elongated plasma by providing $\boldsymbol {B}^\text {coil}_{t_1}$. Regardless of the magnetic field updating scheme, when ${\boldsymbol {B}^{*}}^\text {eddy}_{t_1} = \boldsymbol {B}^\text {eddy}_{t_0}$ is given as the initial guess, the output of VMEC results in an upwardly displaced equilibrium. This outcome arises because, to restore force balance within the fixed input magnetic field distribution, which is a vertically unstable configuration, the equilibrium must shift upwards to obtain the upward electromagnetic force. Eddy currents induced by this displacement exert a downward force that pushes the plasma back. If the magnetic field is simply updated, the added downward electromagnetic force from the eddy currents yields an equilibrium with an even further upward displacement in the subsequent iteration. A specific example of such divergent behaviour is presented in § 3.2. To achieve consistent results, it is essential to apply a magnetic field update strategy that guides the calculations towards convergence.

2.2. Optimisation of adjustable parameter $\alpha$

This subsection describes how to update the magnetic field using an adjustable parameter.

The problem defined in § 2.1 for achieving convergence is organised using mathematical expressions. For simplicity, the magnetic field produced by the eddy currents at $t_1$ will be denoted as $\boldsymbol {B}$. Note that all variables refer to their values at $t=t_1$ unless otherwise specified. Since the vacuum magnetic field map provided to VMEC should be given at discrete points in cylindrical coordinates, the magnetic field vector array $\boldsymbol {B}$ has components in the $(R, \phi, Z)$ directions and a spatial distribution for the discrete $(R, \phi, Z)$ grid. If the number of grid points is $n$, it can be expressed as $\boldsymbol {B} \in \mathbb {R}^{3n}$. The magnetic field vector at a spatial point $i$ is denoted by $\boldsymbol {B}_{(i)}$, and the magnitude of its $x$ component ($x$ can be $R$, $\phi$ or $Z$) is further denoted by $B_{(i)x}$.

The sequence of equilibrium and eddy current calculations, as illustrated in figure 1, is represented by the function $\boldsymbol {f}$. The function $\boldsymbol {f}$ takes $\boldsymbol {B}$ as its input variable and outputs the magnetic field, $\boldsymbol {f}[\boldsymbol {B}]$, produced by the newly calculated eddy currents. Here, $\boldsymbol {f}$ is a mapping from a magnetic field distribution vector of size $3n$ to another magnetic field distribution vector of the same size, $\boldsymbol {f} \colon \mathbb {R}^{3n} \mapsto \mathbb {R}^{3n}$. Note that this mapping is nonlinear and includes MHD equilibrium and eddy current calculations. If the input magnetic field distribution $\boldsymbol {B}$ is equal to the output magnetic field distribution $\boldsymbol {f}[\boldsymbol {B}]$, the whole calculation is regarded as being converged, where the self-consistent solution of the MHD equilibrium and eddy current distribution is obtained. In other words, obtaining a self-consistent solution can be interpreted as a problem of finding an input $\boldsymbol {B}$ such that it is a ‘fixed point’ that is not changed by the mapping $\boldsymbol {f}$. To make the problem numerically solvable, it is transformed into an optimisation problem. We define the variation in the magnetic field distribution by $\boldsymbol {f}$ as

(2.1)\begin{equation} \boldsymbol{g}\left[\boldsymbol{B}\right] \equiv \boldsymbol{f}\left[\boldsymbol{B}\right] - \boldsymbol{B}, \end{equation}

which serves as the residual vector array. The aim is to bring this $\boldsymbol {g}[\boldsymbol {B}]$ closer to the zero vector.

Since $\boldsymbol {f}$ includes equilibrium calculation and exhibits nonlinearity, it is impractical to find an exact solution in a single step. We consider an iterative approach that gradually converges to an exact solution.

Here, we define an adjustable parameter $\alpha$ and propose a simple method to optimise it. Due to the large number of dimensions of the input and output variables, it is difficult to determine the dependence of $\boldsymbol {g}[\boldsymbol {B}]$ on $\boldsymbol {B}$, so the usual gradient method cannot be applied. Therefore, with the provisional magnetic field distribution as $\boldsymbol {B}^*$, the arbitrary magnetic field variation as $\boldsymbol {b}$ and the adjustable parameter as $\alpha$, the magnetic field update shall proceed by finding a value of $\alpha$ that makes $\boldsymbol {g}[\boldsymbol {B}^*+\alpha \boldsymbol {b}]$ smaller. For example, if $\boldsymbol {b}$ is a variation which makes $\boldsymbol {g}$ smaller, then $\alpha$ should take a positive value; conversely, if the variation makes $\boldsymbol {g}$ greater, $\alpha$ should take a value of zero or negative. The specific method for determining $\boldsymbol {b}$ will be discussed in § 2.3, so we will not go into detail here.

The following describes how to estimate $\boldsymbol {g}[\boldsymbol {B}^*+\alpha \boldsymbol {b}]$ based on the idea of the Gâteaux differential (Hazewinkel Reference Hazewinkel1995). Given linear topological spaces $X$ and $Y$ and a map $F \colon X\mapsto Y$, the Gâteaux differential in the direction $h \in X$ at $x_0 \in X$ is expressed as follows:

(2.2)\begin{equation} \delta F\left[x_0; h\right] = \lim_{\epsilon \rightarrow 0} \frac{F\left[x_0 + \epsilon h\right] - F\left[x_0\right]}{\epsilon}. \end{equation}

The parameter $\delta F[x_0; h]$ on this left-hand side is called the Gâteaux variation, and the Gâteaux variation has homogeneity of degree 1 with respect to $h$, namely

(2.3)\begin{equation} \delta F\left[x_0; \alpha h\right] = \alpha \delta F\left[x_0; h\right] \end{equation}

holds for any scalar $\alpha$. If we apply the Gâteaux differential to the current problem, we obtain the following expression:

(2.4)\begin{equation} \delta \boldsymbol{g}\left[\boldsymbol{B}^*; \boldsymbol{b}\right] = \lim_{\epsilon \rightarrow 0} \frac{\boldsymbol{g}\left[\boldsymbol{B}^* + \epsilon\boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right]}{\epsilon}. \end{equation}

Here, the Gâteaux variation $\delta \boldsymbol {g}$ is just a vector of the $\mathbb {R}^{3n}$ directional derivatives of each real-valued components of $\boldsymbol {g}$ along the direction of the specific vector $\boldsymbol {b}$. Now consider approximating the limit operation by a finite-width difference operation, as is done in general numerical differentiation. Using an appropriate positive constant $\epsilon _{c}$, we denote the approximated Gâteaux variation as follows:

(2.5)\begin{equation} \Delta \boldsymbol{g}\left[\boldsymbol{B}^*; \boldsymbol{b}\right] \equiv \frac{\boldsymbol{g}\left[\boldsymbol{B}^* + \epsilon_{c}\boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right]}{\epsilon_{c}}. \end{equation}

This approximate Gâteaux variation (2.5) will reproduce the original Gâteaux variation (2.4) if $\epsilon _{c}$ is small enough. Assuming that $\epsilon _{c}$ is small enough to make $\delta \boldsymbol {g} \approx \Delta \boldsymbol {g}$ hold for $\boldsymbol {b}$ and $\alpha \boldsymbol {b}$, the following approximation holds due to the homogeneity shown in (2.3):

(2.6)\begin{equation} \boldsymbol{g}\left[\boldsymbol{B}^* + \alpha\epsilon_{c}\boldsymbol{b}\right] \approx \boldsymbol{g}\left[\boldsymbol{B}^*\right] + \alpha\left(\boldsymbol{g}\left[\boldsymbol{B}^* + \epsilon_{c}\boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right]\right). \end{equation}

The implication of this expression is that the variation of $\boldsymbol {g}$ is approximately homogeneous for small changes in input. Let us now consider the case where $\epsilon _c=1$ in order to make (2.6) an expression for $\boldsymbol {g}[\boldsymbol {B}^* + \alpha \boldsymbol {b}]$

(2.7)\begin{equation} \boldsymbol{g}\left[\boldsymbol{B}^* + \alpha\boldsymbol{b}\right] \approx \boldsymbol{g}\left[\boldsymbol{B}^*\right] + \alpha\left(\boldsymbol{g}\left[\boldsymbol{B}^* + \boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right]\right). \end{equation}

This condition can be also interpreted as that the exact derivative is replaced by a finite differences approximation. This formula becomes less accurate the larger $\boldsymbol {b}$ is, but as a rough approximation it can be considered to hold for any $\boldsymbol {b}$. Even if $\boldsymbol {b}$ is so large that the approximation is completely invalid, (2.7) can be made valid by redefining $\boldsymbol {b}$ as a reduction to a suitable scale. Gathering terms in (2.7) with respect to $\alpha$, we obtain

(2.8)\begin{gather} \boldsymbol{f}_0 \equiv \boldsymbol{g}\left[\boldsymbol{B}^*\right] = \boldsymbol{f}\left[\boldsymbol{B}^*\right] - \boldsymbol{B}^*, \end{gather}

(2.9)\begin{gather}\boldsymbol{f}_1 \equiv \boldsymbol{g}\left[\boldsymbol{B}^* + \boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right] = \boldsymbol{f}\left[\boldsymbol{B}^* + \boldsymbol{b}\right] - \boldsymbol{f}\left[\boldsymbol{B}^*\right] - \boldsymbol{b}, \end{gather}

(2.10)\begin{gather}\boldsymbol{g}\left[\boldsymbol{B}^* + \alpha\boldsymbol{b}\right] \approx \boldsymbol{f}_0 + \alpha \boldsymbol{f}_1. \end{gather}

Based on (2.7), the value of $\alpha$ will be determined such that $\boldsymbol {g}[\boldsymbol {B}^* + \alpha \boldsymbol {b}]$ becomes smaller. Here, we propose two policies. One policy is to aim for the average value of each component going to zero, and the other is to aim for the mean of the sum of squares of the components being a local minimum.

First, the policy to reduce the average value of each component to zero is considered. If expression (2.10) is separated in components, it can be rewritten as

(2.11)\begin{equation} g_x\left[\boldsymbol{B}^* + \alpha\boldsymbol{b}\right] \approx f_{0,x} + \alpha f_{1,x}. \end{equation}

This expression is simply a reformulation of expression (2.10), without introducing any new assumptions. From (2.11), we can find $\alpha$ such that the average value of each component of the residual becomes zero, but such an $\alpha$ will be different for each component. Since it is necessary to reduce all components to bring $\boldsymbol {g}$ close to zero, it is desirable to optimise each component simultaneously. Therefore, we allow $\alpha$ to take different values for each component and consider treating $\alpha$ as a vector

(2.12)\begin{equation} \boldsymbol{\alpha}_\text{vector} \equiv \sum_{x=R,\phi,Z} \alpha_x \boldsymbol{e}_x, \end{equation}

where $\boldsymbol {e}_x$ is the unit vector for the $x$-direction. When using $\boldsymbol {\alpha }_\text {vector}$, the initially assumed updated magnetic field $\boldsymbol {B}^* + \alpha \boldsymbol {b}$ is modified to $\boldsymbol {B}^* + \boldsymbol {\alpha }_\text {vector}\odot \boldsymbol {b}$, where the symbol $\odot$ denotes component-by-component multiplication. One might question the modification of the starting formula in the discussion, but we have merely relaxed the implicit constraint that $\alpha$ is a scalar. To adapt $\boldsymbol {\alpha }_\text {vector}$, the expression (2.7) can be modified to a component-wise approximation formula

(2.13)\begin{equation} g_x\left[\boldsymbol{B}^* + \alpha_x b_x\boldsymbol{e}_x\right] \approx g_x\left[\boldsymbol{B}^*\right] + \alpha_x\left(g_x\left[\boldsymbol{B}^* + b_x\boldsymbol{e}_x\right] - g_x\left[\boldsymbol{B}^*\right]\right), \end{equation}

where we remind the reader that $x$ stands for $(R, \phi, Z)$. Using the symbol $odot$, it reads

(2.14)\begin{align} \boldsymbol{g}\left[\boldsymbol{B}^* + \boldsymbol{\alpha}_\text{vector}\odot \boldsymbol{b}\right] & \approx \boldsymbol{g}\left[\boldsymbol{B}^*\right] + \boldsymbol{\alpha}_\text{vector}\odot\left(\boldsymbol{g}\left[\boldsymbol{B}^* + \boldsymbol{b}\right] - \boldsymbol{g}\left[\boldsymbol{B}^*\right]\right) \nonumber\\ & = \boldsymbol{f}_0 + \boldsymbol{\alpha}_\text{vector}\odot\boldsymbol{f}_1. \end{align}

To make the spatially averaged residuals for each $x$ component vanishing, $\alpha _x$ is determined by

(2.15)\begin{equation} \alpha_x ={-}\frac{\overline{f_{0,x}}}{\overline{f_{1,x}}}. \end{equation}

Here, the overline denotes averaging over the spatial distribution. If $S_\text {in}$ denotes the region for averaging and $n_\text {in}$ denotes the number of spatial points in $S_\text {in}$, then $\bar {\boldsymbol {B}}\equiv ({1}/{n_\text {in}})\sum _{i\in S_\text {in}}{\boldsymbol {B}_{(i)}}$. Since applying independent operations to each component may result in a magnetic field distribution that is not divergence free, it is necessary to address this issue in the implementation to obtain an appropriate magnetic field distribution. This point will be discussed in § 2.3.

In general, independence of each component of $\boldsymbol {g}$ is not exactly established, but a rough component-by-component correspondence is thought to exist. Since changes in the input magnetic field move the plasma and the eddy currents are induced to suppress it, it can be inferred that there is a correspondence between the changes in input and output of $\boldsymbol {f}$ in each spatial direction. Furthermore, because the cylindrical coordinate system is orthogonal, the direction of a spatial vector can be represented by a combination of independent components. From these considerations, it can be assumed that the solution by (2.15) brings $\boldsymbol {g}$ closer to the zero vector, even though strict independence is not established.

Next, the policy of minimising the sum of squares of the components is considered. The spatial average of the sum of squares of $\boldsymbol {g}$ is denoted $\mathcal {L}^2$ as follows:

(2.16)\begin{equation} \mathcal{L}^2\left[\boldsymbol{B}\right] \equiv \overline{{\lVert \boldsymbol{g} \lVert}^2}\left[\boldsymbol{B}\right] = \frac{1}{n_\text{in}}\sum_{i\in S_\text{in}} \sum_{x=R,\phi,Z} g_{(i)x}^2\left[\boldsymbol{B}\right]. \end{equation}

Then, $\mathcal {L}^2[\boldsymbol {B}^* + \alpha \boldsymbol {b}]$ can be estimated using (2.10) as

(2.17)\begin{equation} \mathcal{L}^2\left[\boldsymbol{B}^* + \alpha \boldsymbol{b}\right] \approx \overline{{\lVert \boldsymbol{f}_0 + \alpha\boldsymbol{f}_1 \lVert}^2} = \overline{{\lVert \boldsymbol{f}_0 \lVert}^2} + \alpha^2\overline{{\lVert \boldsymbol{f}_1 \lVert}^2} + 2\alpha\overline{\boldsymbol{f}_0 \boldsymbol{\cdot} \boldsymbol{f}_1}, \end{equation}

and the derivative for $\alpha$ is

(2.18)\begin{equation} \frac{\partial}{\partial \alpha}\mathcal{L}^2\left[\boldsymbol{B}^* + \alpha \boldsymbol{b}\right] \approx 2\alpha\overline{{\lVert \boldsymbol{f}_1 \lVert}^2} + 2\overline{\boldsymbol{f}_0 \boldsymbol{\cdot} \boldsymbol{f}_1}. \end{equation}

From (2.18), $\mathcal {L}^2$ takes a local minimum value at

(2.19)\begin{equation} \alpha_\text{scalar} ={-}\frac{\overline{\boldsymbol{f}_0 \boldsymbol{\cdot} \boldsymbol{f}_1}}{\overline{{\lVert \boldsymbol{f}_1 \lVert}^2}}, \end{equation}

which is the solution for the second policy. Although this solution is equivalent to what has been described briefly in a reference work (Cianciosa et al. Reference Cianciosa, Isernia, Rubinacci, Terranova and Villone2021), here, we have given its derivation more consistently. There could be a variety of theoretically possible optimisation policies other than the two proposed here, but we did not find any very promising ones as far as we tried.

The described method can be interpreted as a single-variable optimisation approach for the parameter $\alpha$, representing a step size in a specific search direction, using Newton's method. The policy of treating $\alpha$ as a vector constitutes a root-finding problem for the spatially averaged values of each component, whereas treating $\alpha$ as a scalar represents a least-squares problem. This approach may also resemble the secant method, as it employs a secant defined by two points, $\boldsymbol {B}^*(\alpha =0)$ and $\boldsymbol {B}^*+\boldsymbol {b}(\alpha =1)$; however, since the correspondence between the value of $\alpha$ and the actual magnetic field distribution varies with each iteration, the proposed method differs from the conventional secant method, which uses the result from the previous iteration to give one of the points defining the secant taken from the result of the previous iteration.

2.3. Implementation

This subsection explains the specific implementation of the method described in § 2.2 into numerical calculations.

In this study, the variation $\boldsymbol {b}$ in the input magnetic field is given by the residual vector for the provisional magnetic field

(2.20)\begin{equation} \boldsymbol{b} = \boldsymbol{g}\left[\boldsymbol{B}^*\right], \end{equation}

as the search direction in which some variation in $\boldsymbol {g}$ is expected, whether positive or negative. Since $\boldsymbol {g}[\boldsymbol {B}^*]$ is required for determining $\alpha$ in any case, it can be used for both purposes to reduce the number of times $\boldsymbol {f}$ needs to be evaluated. The property of $\boldsymbol {g}[\boldsymbol {B}^*]$ becoming smaller as the convergence progresses is convenient because it automatically scales $\boldsymbol {b}$ appropriately in iterative calculations, and leads to a better approximation accuracy of the approximation (2.7). This approach to determining $\boldsymbol {b}$ does not guarantee the most effective reduction of $\boldsymbol {g}$; however, it is practical as it does not require gradient calculations and allows for the straightforward use of VMEC.

Next, the update formula for $\boldsymbol {B}^*$ is considered. Although $\alpha$ was determined in § 2.2 to make $\boldsymbol {g}[\boldsymbol {B}^*+\alpha \boldsymbol {b}]$ smaller, it remains open to consideration whether $\boldsymbol {g}[\boldsymbol {B}^*+\alpha \boldsymbol {b}]$ should be used as it is as an update formula for $\boldsymbol {B}^*$. Since the final solution is the output of $\boldsymbol {f}$, $\boldsymbol {B}^*$ should be a physical solution which is a possible result of the equilibrium and eddy current calculations. It is not guaranteed that a solution equal to $\boldsymbol {B}^*+\alpha \boldsymbol {b}$ will be given by the results of the series of calculations, even if $\boldsymbol {B}^*$ is a physical solution, because $\boldsymbol {B}^*+\alpha \boldsymbol {b}$ is numerically deformed from $\boldsymbol {B}^*$. For example, $\boldsymbol {B}^*+\boldsymbol {\alpha }_\text {vector}\odot \boldsymbol {b}$ is not a physically appropriate solution because it does not keep solenoidality due to the deformation by $\boldsymbol {\alpha }_\text {vector}$. Therefore, especially when $\alpha$ is treated as a vector, the output of $\boldsymbol {f}$ with $\boldsymbol {B}^*+\boldsymbol {\alpha }_\text {vector}\odot \boldsymbol {b}$ as input is adopted as the provisional magnetic field $\boldsymbol {B}^*$ for the next iteration, namely

(2.21)\begin{equation} \boldsymbol{B}^* \leftarrow \boldsymbol{f}\left[\boldsymbol{B}^* + \boldsymbol{\alpha}_\text{vector} \odot \boldsymbol{b}\right]. \end{equation}

On the other hand, the increase in the number of evaluations of $\boldsymbol {f}$ is not desirable for computational efficiency. Considering that $\boldsymbol {B}^* + \boldsymbol {b} = \boldsymbol {f}[\boldsymbol {B}^*]$ is a physical solution when (2.20) is adopted, $\boldsymbol {B}^*+\alpha _\text {scalar}\boldsymbol {b}$ is a mixture of two physical solutions, a distribution that is possible as a result of eddy current calculation (although the existence of equilibrium calculation results that induce such eddy currents is still not guaranteed). Therefore, when $\alpha$ is a scalar, $\boldsymbol {B}^*+\alpha _\text {scalar}\boldsymbol {b}$ has a small deviation from physical solutions, and the following update formula can be used:

(2.22)\begin{equation} \boldsymbol{B}^* \leftarrow \boldsymbol{B}^* + \alpha_\text{scalar} \boldsymbol{b}. \end{equation}

Comparing formulae (2.21) and (2.22), it should be noted that treating $\alpha$ as a vector requires one extra evaluation of $\boldsymbol {f}$. Figure 2 shows a flowchart based on the implementation methods described so far, which represents the contents of one iteration of the convergence calculation.

Figure 2. Flowchart showing one iteration of the convergence calculation to update the provisional vacuum magnetic field distribution $\boldsymbol {B}^*$.

One point to note in the implementation is the numerical divergence in expression (2.15). On the right-hand side of (2.15), there are cases where the denominator is almost zero but the numerator is a finite value, making $\alpha _x$ extremely large. To avoid such a problem, it is recommended to provide an appropriate upper limit for the value of $\alpha _x$ (e.g. a fixed value or a value set with reference to the value of $\alpha _\text {scalar}$).

3.1. Horizontal displacement

This subsection deals with the case where only horizontal displacements are induced by the reduction of the total toroidal plasma current.

3.1.1. Results for simple update

We first illustrate the results obtained when using a simple update method, where the output magnetic field from the previous iteration is directly used as the input for the next iteration, expressed by the substitution equation as $\boldsymbol {B}^* \leftarrow \boldsymbol {f}[\boldsymbol {B}^*]$, without employing the method described in § 2.

Figure 5 shows the results of iterative calculations for 10 iterations. In figure 5(b), which shows the change in the horizontal position of the magnetic axis, the results using the proposed method in § 2 is also plotted for comparison (these results will be discussed in § 3.1.2). One iteration means a series of calculations to update the vacuum field distribution. The label ‘Simple update’ in figure 5(b) represents the results obtained using the simple magnetic field update method, corresponding to the outcome shown in figure 5(a). The labels ‘$\boldsymbol {\alpha }_\text {vector}$’ and ‘$\alpha _\text {scalar}$’ refer to the series of calculations shown in figure 2.

Figure 5. Changes in the plasma shape and the position of magnetic axis during a simple iterative calculation with alternating updates of equilibrium and eddy currents, following a 1 % decrease in total toroidal current. (a) Shows the plasma surfaces. The grey line represents the conductor simulating the VV with flowing eddy currents. The dashed line shows the initial equilibrium plasma shape before the current decreases, while the solid coloured lines show the equilibrium plasma shapes after the current decreases. (b) Shows the horizontal positions of the magnetic axis during the iterations. The green line is the result of simple updates, and the blue and orange lines are examples of the result using the proposed method.

From figure 5(a), it can be seen that the plasma is significantly displaced and deformed by the iterative calculations with simple update, and the calculation results do not converge. ‘Simple update’ in figure 5(b) shows that the amplitude is increasing, with the direction of displacement reversing in each iteration. Focusing on the 0th and 1st iterations in 5(b), at the 0th iteration, all results show the inward displacement due to the current decrease, but at the 1st iteration, the magnetic axis of ‘Simple update’ is located obviously outside the initial position, while ‘$\boldsymbol {\alpha }_\text {vector}$’ and ‘$\alpha _\text {scalar}$’ quickly settle near the final convergence position. Based on the equilibrium displaced outward from the initial position (the position at the previous time step), the eddy currents are induced to push the plasma back inward. As a result, in the 2nd iteration, adopting the magnetic field produced by these eddy currents, the plasma becomes even more excessively displaced inward compared with the 0th iteration. In other words, during the convergent iteration of evolution involving horizontal displacement, although the initial displacement occurred in a physically reasonable direction, the excessive magnitude of the displacement led to oscillatory divergence in the behaviour.

The results in figure 5 demonstrate that simple updates cannot correctly calculate the evolution of the equilibrium considering eddy currents, although the simulation set-up is very simple. Since the same behaviour was observed even when the time step or the amplitude of current decay was changed, this behaviour is considered to be a qualitative instability caused by the computational scheme rather than numerical instabilities (e.g. instability due to the time step width). Given that this behaviour subsides when the resistivity of the VV is increased or the distance between the plasma and VV is enlarged, such behaviour should be specific to situations where a sufficiently conductive wall is present.

3.1.2. Results for the proposed method

The results of using the proposed method in the current reduction simulation are presented here and the differences due to the treatment of $\alpha$ are discussed.

Figure 6 visualises the histories of the residual metric $\mathcal {L}^2[\boldsymbol {B}^*]$ defined in (2.16) during the convergence calculation. In both ‘$\boldsymbol {\alpha }_\text {vector}$’ and ‘$\alpha _\text {scalar}$’, $\mathcal {L}^2$ gets smaller with iterations, reaching a value of order $1\times 10^{-20}$, which appears to be due to numerical errors. Although there appears to be no significant difference in the global convergence speed, the actual computational cost is smaller for ‘$\alpha _\text {scalar}$’ because one extra evaluation of $\boldsymbol {f}$ is required in ‘$\boldsymbol {\alpha }_\text {vector}$’, as shown in figure 2.

Figure 6. History of decrease in $\mathcal {L}^2$ of the magnetic field during convergence calculation for temporal development with the plasma toroidal current decrease. The label $\boldsymbol {\alpha }_\text {vector}$ represents that $\alpha$ is treated as a vector and $\alpha _\text {scalar}$ represents that the parameter $\alpha$ is treated as a scalar.

3.2. Combination of horizontal and vertical displacements

This subsection examines the calculations in the presence of vertical displacement in addition to the horizontal one.

By shifting the VV by ${+0.1}\ {\rm m}$ in the $Z$-direction, the relative initial position of the plasma is shifted downwards. This operation can be interpreted as placing the initial equilibrium plasma below the neutral equilibrium location for the vertical displacement events (Nakamura et al. Reference Nakamura, Yoshino, Neyatani, Tsunematsu, Azumi, Pomphrey and Jardin1996). The physical mechanism is that the eddy currents induced by the reduced plasma current are in the same direction as the plasma current, so the attraction between the same direction currents should cause the plasma to displace downwards (nearer side to VV), while also inducing eddy currents to suppress it. Naturally, there is also the inward shift of the plasma due to current damping, as in § 3.1, so the situation is complicated by a combination of horizontal and vertical displacements.

3.2.1. Results for simple update

As in § 3.1.1, the result of the case without using the proposed method is checked first.

In this case, blend update is applied instead of the simple update to suppress the horizontal displacement. Otherwise, the calculation would collapse due to increasing horizontal displacements as in figure 5 before the vertical displacements become significant. In blend update, the current provisional magnetic field is blended to the new magnetic field when updating the magnetic field as $\boldsymbol {B}^* \leftarrow (\boldsymbol {B}^* + \boldsymbol {f}[\boldsymbol {B}^*])/2$. In order to ensure that a visible change appears, in this calculation the total toroidal plasma current decrease was set to 5 % instead of 1 %.

Figure 7 shows the history of the plasma shape and the vertical position of the magnetic axis in convergence calculations. From figures 7(a) and 7(b), it can be seen that the vertical displacement increases monotonically for the blend update without using the proposed method.

Figure 7. Changes in the plasma shape and the position of the magnetic axis during a simple iterative calculation with alternating updates of equilibrium and eddy currents, adopting blend update, following a 5 % decrease in total toroidal current. The VV is shifted ${+0.1}\ {\rm m}$ in the $Z$-direction from the position shown in figure 4. (a) Shows the plasma surfaces. (b) Shows the vertical positions of the magnetic axis during the iterations.

Focusing on the direction of displacement, the displacement is occurring in the upward direction, contrary to the physical mechanism. This behaviour is thought to be exactly due to the mechanism described in § 2.1. It is shown that, by repeating equilibrium calculations with a fixed input magnetic field distribution and eddy current calculations based on the equilibrium plasma, displacements in the direction opposite to the electromagnetic force are amplified.

What the results in § 3.1.1 and this section show is that the relationship between plasma displacement and eddy currents can vary with direction, in the VMEC calculation. In equilibrium calculations, the plasma shape is determined such that the residual forces acting on the plasma are minimised. However, in calculations where the input magnetic field is fixed, the displacement direction that reduces the residual force corresponds to the same direction as the residual force in the stable direction and the opposite direction in the unstable direction. Physically, eddy currents are induced to exert forces that suppress the displacement from the previous time step, and the magnetic field produced by the eddy currents works to suppress displacement in the stable direction, as shown in figure 5, and amplify it in the unstable direction in the next equilibrium calculation, as shown in figure 7.

3.2.2. Results for the proposed method

The results of the proposed method are discussed in this section. To keep consistency with the results in figure 7, the current decrease is set to 5 % as in § 3.2.1.

Figure 8 visualises the histories of $\mathcal {L}^2[\boldsymbol {B}^*]$ during the convergence calculation. Whereas $\mathcal {L}^2$ decreases with each iteration and reaches the order of numerical error as in figure 6 for ‘$\boldsymbol {\alpha }_\text {vector}$’, the decrease in $\mathcal {L}^2$ almost stops progressing at of the order of $1\times 10^{-8}$ for ‘$\alpha _\text {scalar}$’. As can be seen in figure 7, the vertical position of the magnetic axis in ‘$\alpha _\text {scalar}$’ oscillated around zero, and the vertical displacement did not seem to converge to the correct result. This behaviour can be explained by the different characteristics depending on the direction of displacement, also mentioned in § 3.2.1. In the horizontal direction, the eddy currents act in the direction of suppressing displacement. The residual vector, which is the change from input to output of $\boldsymbol {f}$, is therefore a magnetic field in the direction of suppressing the displacement of the plasma, and the sign of the appropriate $\alpha _\text {scalar}$ in formula (2.22) is usually positive. In the vertical direction, on the other hand, the eddy currents act in the direction of amplifying the displacement, so the appropriate $\alpha _\text {scalar}$ is negative. In ‘$\alpha _\text {scalar}$’, the magnetic field is updated using a single $\alpha _\text {scalar}$ regardless of the direction, so an update that simultaneously suppresses displacements in both directions is basically impossible. In particular, ‘$\alpha _\text {scalar}$’ is based on local minimisation of $\mathcal {L}^2$, so it tends to stay in the local solution, and no matter how far the iterations are increased, it is unlikely to move towards the correct solution.

Figure 8. History of decrease in $\mathcal {L}^2$ of the magnetic field during convergence calculation for temporal development with the toroidal plasma current decreasing.

From this result, it can be concluded that the policy of treating $\alpha$ as a scalar may not be effective and that it is preferable to adopt the policy of treating $\alpha$ as a vector. inputsections/application