Expanded capabilities of the CarMa code in modeling resistive wall mode dynamics with 3-D conductors

In this work, an improved version of the CarMa code is presented, called CarMa-D, for the analysis of resistive wall modes (RWMs) in fusion devices, simultaneously considering the effects of volumetric three-dimensional conducting structures, plasma dynamics, toroidal rotation or drift-kinetic damping. Unlike static CarMa, CarMa-D does not rely on the simplifying assumptions such as neglecting the plasma mass, toroidal rotation and kinetic damping. The new coupling strategy is based on matrix-based Padé rational functions approximation of plasma a response. The arising mathematical model is formally equivalent to the original CarMa model, but with a higher number of degrees of freedom to model the dynamics of the plasma. CarMa-D overcomes the main limitations of the original CarMa, in particular: (i) the massless assumption for the plasma is removed, allowing modeling of global modes growing both on ideal kink time scales and in the typical RWM growth rate regime, with a suitable treatment of the model; (ii) the effects of toroidal plasma flow and drift kinetic damping can be included into the new model, providing a powerful tool to study macroscopic phenomena where both plasma dynamics and 3-D conducting structures play important roles.


Introduction
One of the most stringent restrictions to reach high plasma performances in fusion devices is given by external ideal MHD modes of low toroidal mode number n [1,2], driven either by the plasma current or pressure. These instabilities produce magnetic field perturbations, which induce stabilizing image currents in the metallic structures placed around the plasma. Such image currents, flowing inside the resistive wall, reduce the mode growth rate from the Alfvén time scale (microseconds) to the characteristic time of the passive conducting structures (the so called wall time, i.e. the penetration time of the perturbed magnetic field through the wall). This time ranges from several milliseconds for the existing devices to about 0.3 s for ITER [3,4]. The residual kink instability is called resistive wall mode (RWM), which grows at much slower time scale with respect to the original ideal kink, opening to the possibility of stabilization with magnetic feedback.
It is well known that the three-dimensional characteristics of conducting structures, as well as the geometry of feedback coil system, play an important role in determining the mode dynamics [5]. It is therefore important to implement modelling tools that take these characteristics into account. Many different strategies have been applied to tackle active control of RWMs. Both experimental and modelling examples can be found on a variety of Tokamak devices [6][7][8][9], and extensions of the virtual shell concept for stabilization of multiple RWMs have been applied to reversed field pinch experiments as well [10][11][12][13].
The CarMa code [14,15] is a computational tool developed for the purpose of studying RWM instabilities coupled with conducting structures by self-consistently interfacing the CAR-IDDI code [16,17] with the MARS-F code [18]. CarMa takes rigorously into account the real geometry of conducting structures, such as the thickness of the conductors. Many other tools have been developed with similar purpose [19][20][21][22][23], but mainly using a thin wall approximation [24] for the passive conductors. The CarMa code overcomes this assumption, taking advantage of a volume integral solution of the eddy current problem provided by CARIDDI. This way, it is possible to evaluate the effect on RWMs of conductors in which the penetration depth is comparable with their width [25].
CarMa has been successfully used in existing devices, such as JET [26] and RFX-mod [5], with applications to RWM control modelling [27] and control system optimization [28]. Predictive simulations have been carried out for future devices as well, such as JT-60SA [29] and ITER [25,30,31]. However, the coupling strategy behind the CarMa code suffers of certain limitations: the main one is the assumption of static plasma response to external magnetic field perturbation, the so-called massless approximation. This means that the plasma inertia and any associated Alfvén-wave-like phenomena are neglected on the time scale of interest [32].
This work reports an essential improvement that we have made to CarMa. This consists in abandoning the massless approximation and implementing a proper treatment of the plasma dynamic effect, by introducing matrix based rational function approximation for the plasma response. The new tool allows self-consistent coupling of magnetic feedback with plasma dynamic effects such as the plasma toroidal rotation and the driftkinetic damping on the RWM. This latter aspect in particular allows coupling with the MARS-K code [33]. The improved version of the code is named CarMa-Dynamic, or CarMa-D. The paper is organized as follows. Section 2.1 briefly summarizes the original CarMa code, reporting the key equations of the coupling strategy. This is followed by section 2.2 describing CarMa-D, where details of the matrix based rational interpolation of the plasma response are reported and the resulting mathematical model is formulated. In section 3, the new code is tested against a n=1 kink instability for a circular shaped plasma. Conclusions and final remarks are reported in section 4.

Standard CarMa coupling based on static plasma response
The CarMa code is based on the control surface concept to self-consistently couple a linear MHD solver (MARS-F [18]) to a 3-D code for the computation of eddy currents in the metallic structures surrounding the plasma (CARIDDI [16,17]): a solid theoretical formulation of the coupling procedure can be found in [14].
The equilibrium is assumed to be toroidally symmetric: the region of space occupied by the plasma in its reference equilibrium is called V p and is enclosed by a suitable coupling surface S e , which does not intersect neither the plasma boundary nor the conducting structures, and does not enclose any conducting structures inside the surface (see figure 1 for the reference geometry). This coupling surface is the key feature of CarMa coupling strategy: the electromagnetic effects of the plasma, seen from the external environment, is described as produced by an equivalent current density j eq flowing on the coupling surface [34].
The most important assumption at the base of the aforementioned coupling procedure is that the plasma is assumed to be static, i.e. the plasma mass is neglected (massless approximation) and no toroidal flow is considered (plasma equilibrium fluid velocity = V 0 0 ). Under this hypothesis, the governing equations inside the plasma volume V p are single fluid, linearized MHD equations: where P is the plasma pressure, B  is the magnetic flux density, J  is the current density, x  is the plasma displacement and Γ is the specific heat ratio. The subscript 0 is for equilibrium quantities, while 1 means the first order variation. With respect to the usual formulation of the single-fluid ideal MHD equations, the terms related to plasma mass have been neglected. The equations for outside vacuum region V e bounded by the surface S e are: At this point, the components of the perturbed magnetic field normal and tangential to the coupling surface S e are expanded in poloidal and toroidal Fourier harmonics. For the mode analysis, only one specific toroidal harmonic component (e.g. n = 1) and a suitable spectrum of M poloidal harmonics is considered. It is clear that, in principle, an infinite number M of poloidal harmonics is required: for a realistic implementation only a finite number M is chosen, high enough to ensure the desired accuracy for a given shaped equilibrium. The number of harmonics has been found performing a convergence study on the eigenvalue. In the next part of the work, the lowercase letter (i.e. b) will be used for quantities defined in Fourier space.
The key steps developed in [14], and here briefly reported, give the relations needed to compute the plasma and external contribution, respectively b N pla and b N ex , to the total perturbed magnetic field b N , that are normal to the coupling surface S e . In particular, the plasma contribution is given in terms of equivalent surface current density j eq on the coupling surface instead of b N pla : where the tilde means a vacuum quantity (i.e. without plasma). This mean that K is the matrix mapping the normal field b N imposed on S e to the tangential field b T produced by the plasma, and K is obtained imposing the same b N on S e , but solving vacuum equations (6)- (7). From their combination the matrix F can be obtained, which gives, for a certain perturbation b N on S e , the response of the plasma in terms of equivalent surface current density. These relations have dimension M×M. The second key relation involves b N ex , the external contribution (due to the passive conductors and the feedback actuators) to the total perturbed magnetic field: where W is the M×M static relation which provides b N ex starting from the total field b N . Since we are interested in the study of the RWMs stabilization, the external perturbation b N ex comes from both excitation coils and the eddy current induced inside the conducting structure by the mode itself. The computational problem in the conductors outside the coupling surface is treated through an integral formulation [16], which assumes as primary unknown the current density in conducting structures. The outcoming equation for the eddy current dynamics is: where i represents the eddy currents, R and L are the resistance and self-inductance matrices of the passive conductors, respectively, D is an incidence matrix between the nodes where the voltage is applied and the edges mesh elements, u is the magnetic flux produced by the plasma and linked with the passive conductors, thereforeu d dt is the voltage induced by the plasma evolution and v is the vector of voltages applied to the feedback coils. The suitable manipulation of (8) and (9) allows to write the plasma contribution u d dt in equation (10) in terms of the so called plasma response matrix X, which can be considered as a modification of the inductance operator L: here L* is the perturbed inductance operator which takes into account also the plasma contribution in the eddy current equation.
It is worth noting that X can be defined only with the assumption of inertia-free plasma, the first important assumption behind this coupling scheme: this is reasonable if we are interested in analyzing phenomena which are orders of magnitude slower than Alvén timescale. This is the case of the typical RWM regime, in which the growth/damping rate has the same order of magnitude as the wall characteristic time, a common quantity for the study of RWMs stability [35].
The other assumption is that a fluid description is adopted, i.e. no kinetic damping physics is considered. If the damping physics resulting from the plasma flow and the drift kinetic resonance effects has to be considered, or, conversely, if the dynamical effects play a crucial role, the massless approximation is no longer suitable. In the following, some specific cases will be discussed in which the massless approximation is not valid and hence a suitable extension of the procedure is needed.

New
CarMa-D coupling based on dynamic plasma response 2.2.1. Dynamic plasma response and matrix-based rational function approximation. The inertia-free approximation, from the point of view of the matrix K already introduced in Subsection 2.1, leads to a static relation between the normal and tangential components of the perturbed magnetic field on the coupling surface. When the plasma is close to the idealwall beta limit, or when the plasma flow effects are important, such assumption is no longer valid, and all the quantities involved in the plasma response matrix computation in section 2.1 should depend on the dynamics of the perturbation field, hence on the plasma dynamics. Thus, working in the Laplace domain, system (1)-(5) can be written as [36][37][38]: where s is the complex Laplace variable, corrected by the Doppler shift W in , Ω is the plasma rotation frequency along the toroidal direction f , R the plasma major radius, ρ is the unperturbed plasma density, and ẑ the unit vector in the vertical direction. Other quantities have already been defined in section 2.1. The set of equations (12)- (16) are the full set of linearized MHD equations, written in term of perturbed velocity V, in which the plasma inertia is not neglected. Drift-kinetic effects are included in the MHD model through the equation involving the perturbed kinetic pressure tensors [33]: where p is the pressure tensor, p is the scalar pressure perturbation, p  and p ⊥ are respectively the components of the kinetic pressure perturbations parallel and perpendicular to the equilibrium magnetic field, I is the unit tensor and = b B B | |. The full pressure tensor p is self-consistently included into the MHD formulation by replacing the term p the momentum equation (12).
From the aforementioned considerations it follows that, if the plasma mass is taken into account, a certain dependance of the matrix K on the complex Laplace variable s is expected, giving rise to a dynamic mapping K s ( ). From a formal point of view, the matrix K s ( ), which has now to be considered a matrix function, can be computed as done in section 2.1. To this purpose, the system of equations (12)-(16) written in a reformulated way [36] and solved for frequency dependent boundary conditions b s N ( ) on the coupling surface: where x is the vector of unknowns, and the matrices U, A are obtained through numerical discretization of the differential equations (12)-(16) [39]. This means that equations (8) and (9) become: where W F s s , ( ) ( ) are now matrix-valued transfer functions that describe, together, the plasma response to externally applied magnetic field perturbation.
The arbitrary dependence of the plasma response with respect to the perturbation frequency can be treated as a dynamical linear system: equations (8), (9) and (10) are written showing the explicit dependence of the response matrices with the respect of s. Assuming vanishing initial conditions we obtain: where the matrix M is the mutual inductance between the passive conductors and the plasma equivalent surface current, and the matrix Q maps the eddy currents into the perturbed magnetic field normal to the coupling surface S e . Every entry of the matrices F s ( ) end W s ( ) is a scalar function of the variable s. This system of equations describes the coupled problem of the plasma response to external field perturbation together with the eddy current equation, without any assumption of neglecting the plasma mass. It would be desirable to have an analytic relation for the matrix function W s ( ), F s ( ). However, in this work, the problem is addressed in a numerical way, looking for a suitable set of matrix-based interpolating functions. For this purpose, we recall a result already proved in [40] for cylindrical geometry, stating that the plasma response is a rational function, also called Padé interpolation, of the complex variable s. The Padé approximation, which has been widely adopted in many fields as a model reduction technique, has been used also for representing the plasma response for feedback stabilization of RWMs [6,36].
The Padé approximant for F W s s , ( ) ( ), are matrix-based functions of s, defined as: where E is the identity matrix and A B C D , , , Since the relations required to obtain the matrix coefficients are the same for (22) and (23) With this approach, the plasma response is approximated through a matrix based rational interpolation of the reference one, computed starting from MHD equations (12)- (16). Dealing with an interpolation problem means that the coefficients A i , B i have to be computed starting from the knowledge of the matrix function P s ( ), which has to be interpolated on a certain number of frequencies s i (the socalled basis points). For a given degree k, the coefficients of (23) require the choice of + k 2 1 basis points. These points consist in a suitable set of complex frequencies chosen along the complex plane, and used as input frequencies of the forcing problem (18): this means that, given i as complex excitation frequency, for each s i the plasma response P s i ( ) has to be computed solving the MHD kinetic-hybrid equations (12)-(16) running MARS-K. This is one of the main differences with the static CarMa code, in which the computation of the (static) plasma response has to be computed just once for vanishing excitation frequency  s 0 | | . Among all possible sets of basis points, a desirable choice can be to force the interpolant to match the reference function P s ( ) for two particular frequencies: , to match the static response with this particular choice, the coefficient A 0 is the same matrix for the original CarMa code; , to match the response at infinite frequency: which can be solved by inverting the matrix Z. It is worth noting that the matrix Z can be very ill-conditioned, or even singular to machine working precision. This rank deficiency can happen if the information in the pairs P s s , i i ( ( )) are linearly dependent, and this can be due to several reasons. One possible case is that the plasma response is almost constant for all the range of frequencies, and this is the case where the massless approximation holds and the static CarMa coupling strategy should be properly used. Another possible issue can be a non appropriate choice of the basis points, that could be erroneously placed in a 'flat region' of the plasma response giving rise to a set of linearly dependent pairs P s s , From a purely formal point of view, the desired interpolation degree k can be arbitrarily high. However, the work described in [40] shows analytically, for a cylindrical plasma, that the system exhibits second order dynamics. In addition to this, the higher the interpolation order k is, the worse is the conditioning of the matrix Z, which can compromise the computation of the coefficients A B , i i . This happens because, when the interpolation order k is too high in comparison with the system dynamics, there is not enough information on the data from the MHD computations to obtain the coefficients.
or equivalently, with a compact formalism: ] is the augmented state vector and = u Dv 0 T [ ] is the input quantities. It is well known that a system of differential equations with degree greater than one can be written With obvious redefinitions, the new systems becomes: which is formally equivalent to equation (15)a in [14], that is the original CarMa code, but with a higher number of states to take into account the plasma dynamics. In particular, both R L , * * are sparse matrices, and their dimensions are +ḱ M N D ) , with N 3D the number of discrete unknowns of the eddy current problem.
As already discussed, if the previous approach is used to study a Resistive Wall Mode in the so-called typical RWM regime, e.g. if the instability growth time has the same order of magnitude of wall time, the matrix L* in equation (30) can be singular to the machine working precision. In the static limit, when the growth rate is small if compared to the wall time, the matrices A B , i i in equation (28), with > i 1, are all zeros. This means that the matrix L* is rank deficient, and the eigenvalues of the system - , which are the growth/ damping rates of the RWM, cannot be computed anymore because L* is not invertible. This issue can be overcome exploiting a suitable treatment of the problem. To this purpose we consider that the matrix R* is always invertible, because it is the composition of the identity matrix with R a , which always has a full rank. A possible solution to this problem is to compute the eigenvalues of the system - , to obtain the growth/damping time of the instability rather than the growth/damping rate: the growth rate then follows as the inverse of the growth time. This is one of the strengths of the proposed approach, that is allowing the treatment of global modes both on ideal kink time scales and in the typical RWM growth rate regime with a general strategy, regardless the possible presence of singular matrices.
It can be proved that an analytical relation is available for the block inversion of - . This latter step is not computationally expensive, because the matrix R a is sparse, as can be seen from equation (28). It follows that:  (32), which is the sub-matrix related to the static CarMa system. This means that the number of non-zero elements of the system matrix is almost due to the passive block, and therefore, even if the number of unknowns is much higher, the computational effort to find the eigenvalues does not grow substantially if the problem is solved through a suitable sparse linear algebra toolbox.
With this formalism, in case of static response, the resulting eigenvalues of (32) would be all zeros, corresponding to the instantaneous static response, except for the N 3D eigenvalues of the static CarMa matrix just mentioned, therefore giving exactly the inverse of the eigenvalues of the static CarMa system.

Numerical verification of models
The CarMa-D model has been tested on a reference case in order to show its validity and effectiveness. The benchmark case is a stability analysis of a tokamak plasma with circular cross section (figure 2(a)), with major radius of R 0 =2 m, minor radius of a=0.4 m. The equilibrium pressure and safety factor profiles are the same of [14]. The coupling surface and the resistive wall shape are also plotted in figure 2(a). A circular resistive wall with thickness δ=5mm placed at a w =1.3a. The penetration time of the n=1 magnetic field perturbation is τ w =4.5 ms. Since the primary goal of this work is the validation of CarMa-D with MARS-K, which uses the thin-wall approximation, only one mesh element along the shell thickness is considered.
The validation of CarMa-D has been performed for two different damping physics mechanisms. The first is the parallel sound wave damping model (SD), described in detail in [41,42], which uses a numerically adjustable coefficient to tune the strength of the damping effect. The second damping   model is the self-consistent, full toroidal drift kinetic-model (KD), presented and validated in [33]. Figure 2(b)-(c) reports the MARS-F/K computed RWM eigenvalues assuming the fluid and drift-kinetic models respectively, while scanning the plasma rotation frequency w E . The latter is assumed to be a constant along the plasma minor radius. CarMa-D computations will be benchmarked against these results.
For both damping mechanisms, validation of CarMa-D follows the same procedure as in [43], i.e. performing a scan of the resistivity η of the vessel. In the limit of infinite resistivity, the no-wall growth rate of the ideal kink mode must be recovered, which has the same order of magnitude of Alvénic timescale, and therefore strongly depends on the dynamical effects. The next subsections show the results of this analysis for both damping mechanisms.

RWM with sound wave damping (SD)
As further proof that the plasma response to external magnetic field perturbation should depend on the frequency of the perturbation, figure 3 shows the matrix function P s ( ), computed for different values of s. As was already said, P s ( ) gives, for a certain perturbation b N on S e , the response of the plasma in terms of equivalent surface current density. In this case, the matrix has been computed using, as boundary conditions, the time-varying normal field b s N ( ). For the analyzed case, the complex frequency, normalized by the Alfvén frequency w A , varies in the range w Î´s 1 10 , 1 10 In figure 3 only some frequencies have been shown. As can be seen, the matricial function is clearly not constant, thus for this case the assumption of static response is not appropriate.
As already said, the function P s ( ) has to be interpolated, thus the reference P s ( ) has to be evaluated on a certain number of frequencies s i through equations (12)- (16). The choice of these frequencies, the so-called basis points, can be done in several ways or along different directions on the complex plane. To test different configurations, we have considered three different paths along the complex plane where the basis points are placed, named respectively complex (s=σ+iω, imposing σ=ω for the computations), real (s=σ) and imaginary (s=iω) cases.
In order to test the effectiveness of rational interpolation, the first important step is to understand which set of basis point is the most suitable to be chosen as complex interpolation frequencies. To do so we can compare the reference values of P s ( ) with its interpolant. The function P s ( ) is a matricial function of s, meaning that every entry p s m m is a scalar function of s. On the other hand, for a chosen s i =γ i +iω i , the quantity P s i ( ) is a M×M matrix. Figure 4 shows the results for 3rd order interpolant of P s ( ) (only the two major entries of the matrix function are shown) for the three different set of basis points. This benchmark results show that the real set of basis points is much more effective than the other sets of points. This is due to the fact that the real part of the unstable eigenvalue is several times larger than the imaginary part, as can be seen in the rotation scan results reported in figure 2(b). Choosing to interpolate the plasma response along the real axis allows to obtain more information then the other choices of basis points. A comparison of the reference P s ( ) with respect to the interpolated one is important to understand if the set of basis points is suitable, and if the coefficients A B , i i have been computed correctly. As already said, to assess the capabilities of CarMa-D, a scan of the resistivity of the conducting wall is performed. Since the real data gives the best results, the detailed results of the resistivity scan are reported only for this case, in figure 5. These plots compare the CarMa-D results of the different interpolation approaches presented in section 2.2.1, more precisely interpolations of 2nd order for F W s s ( ) ( ) and 3rd order for P s ( ). Results from st 1 order interpolation on F W s s ( ) ( ) are also reported, being the first version of the improved CarMa code presented in [43]. Selected results from figure 5 are also reported in table 1, in order to give a quantitative comparison of the CarMa-D accuracy. Figure 5(c) shows the relative error of the CarMa-D results compared to the reference MARS-F results. The relative error is defined as: Results clearly show that interpolating directly the function P s ( ) gives much more accurate results that interpolating F W s s , ( ) ( ) separately, but a higher order, i.e. 3rd order for this equilibrium, is required. Specifically, the 3rd order interpolation on P s ( ) gives an error around 1% for all the range of resistivity.
A straightforward way of further improving the accuracy would be to increase the order of interpolation. This should be true from a theoretical point of view, but there are numerical issues that gives an upper bound to the interpolation order. Indeed, increasing k in equation (23), the condition number of the matrix Z in equation (26) would rise simultaneously, eventually becoming too high to give reliable coefficients for the interpolation. Further comparison had shown that the accuracy is worse if the interpolation is increased from 3rd to 4th order, showing that the 3rd order interpolation of P s ( ) is the best compromise between proper dynamic description and numerical accuracy. For the same considerations, the best results for F W s s , ( ) ( ) are obtained with 2nd order interpolation. To support the analysis of the computational effort in treating the matrix (32), figure 6 shows its sparsity pattern for the k=3 interpolation degree on P s ( ). The matrix dimensions are 5988×5988, and the non-zeros are »4 10 6 , while for standard CarMa the non-zeros are »3.9 10 6 .

RWM with drift kinetic damping (KD)
The same validation procedure reported in figure 4 has been repeated also for the drift kinetic model, showing exactly the same outcomes. As for the previous test case, the real set of basis point again turns out to be the most effective choice. For this reason, the scan of wall resistivity has been performed only for the real set of basis points, and can be seen in figures 7(a)-(c). On the other hand, the effectiveness of interpolating directly P s ( ) rather than F W s s , ( ) ( ) separately seems no longer clear for the this benchmark case: figure 7(c) shows that both strategies give very good results. Figure 8 compares the eigenvectors, i.e. the spatial pattern of the eddy currents which would be induced in the conducting structures by the mode evolution, with both poloidal and toroidal components computed by MARS-K and CarMa-D, assuming 3rd order interpolation on P s ( ) and real set of basis points. Excellent agreement of the mode eigenvectors is obtained. It is worth noting that the circular cross section and the relatively simple shapes of the current eigenfunction are not representative of modern fusion devices (e.g. ITER), but are used here are simple examples; of course, the capabilities of the code would allow the treatment of more realistic geometries. Table 2 compares the eigenvalues for the RWM case, with η/η ref =1 (left), and the ideal kink case, with η/η ref =10 5 (right). The accuracy both for the RWM and for the ideal kink is excellent.
The results reported in figures 7-8 and table 2 are obtained assuming relatively fast plasma flow (ω A /ω A =10 −2 ). We also report results for a case with slow plasma flow. As shown in figure 2(c), the kinetic effects (thermal particle precessional drift resonances with the mode in this case) occur strongly at very slow rotation frequency [44]. Without going through a detailed investigation of various kinetic damping physics, which is beyond the aim of the present work, it is interesting to test the effectiveness of CarMa and CarMa-D for the particular case, with ω E /ω A =6×10 −4 for two different values of wall resistivity. Table 3 shows that static CarMa erroneously predicts a stable mode, while CarMa-D is very accurate for both resistivity      mesh already used for the previous cases. Figure 9 reports results with the same rotation scan as in figure 2(c) for the kinetic case, computed with CarMa-D both for the axisymmetric mesh and for the 3-D mesh, as well as with static CarMa. With the axisymmetric mesh (red), results given by CarMa-D agree well with the reference MARS-K data (black). The presence of holes in the wall leads, in this specific rotation range, to more unstable RWM [45]. It is interesting to notice that CarMa-D (blue) give less pessimistic results than static CarMa (green). Moreover, for a very low toroidal rotation ω E /ω A =5×10 −4 , the mode growth rate predicted by CarMa-D with the 3-D mesh is the same as the axisymmetric case, because, being the mode growth rate very small, the amplitude of the eddy current is also small, giving a weak effect on the mode. With a negligible stabilizing effect, clearly the geometry of the shell does not influence the mode evolution. On the other hand, static CarMa erroneously predicts a stable mode. The imaginary part of the mode eigenvalue is not strongly affected by the 3-D features of the geometry. An example of the eddy current pattern corresponding to the unstable eigenvector for the three-dimensional mesh is shown in figure 10, where the eddy pattern modification due to the holes can be seen.

Conclusion and discussion
CarMa-D, an improved version of the CarMa code, has been presented. This upgraded code is capable of self-consistently taking into account, in addition to the three dimensional conducting structures, also plasma mass, toroidal rotation and kinetic damping physics for the study of Resistive Wall Modes stability. The CarMa-D code can thus also be used to study ideal kink instabilities. Numerical tests show that, independent of the damping mechanisms assumed for the RWM, the new dynamic model recovers well the reference results provided by MARS-F/K. Furthermore, CarMa-D also predicts growth rates for the ideal kink instability, for which the plasma inertia plays an essential role. An example application of CarMa-D to a circular equilibrium with 3-D conducting structures shows that the inclusion of the drift kinetic damping effects in the model obtained by CarMa-D gives less pessimistic results than what gives the static CarMa. In the future, the CarMa-D code will be applied to several cases of interest, such as existing (i.e. ASDEX-U and MAST-U) and future (ITER and JT-60SA) devices. In particular the code allows to simultaneously take into account passive (e.g. kinetic damping) and active (e.g. feedback) RWM stabilizing contributions, with a detailed description of 3-D geometries of active and passive conductors surrounding the plasma.