Implementation of matrix compression in the coupling of JOREK to realistic 3D conducting wall structures
Federico Cipolletta Nina Schwarz Matthias Hoelzl Salvatore Ventre Nicola Isernia Guglielmo Rubinacci Alejandro Soba Mervi Johanna Mantsinen the JOREK Team


Matrix compression of response matrices in JOREK]Implementation of matrix compression in the coupling of JOREK to realistic 3D conducting wall structures





^1 Barcelona Supercomputing Center (BSC), Barcelona, Spain

^2 Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching b. M., Germany
^3 Dipartimento di Ingegneria Elettrica e dell’Informazione, Università degli Studi di Cassino e del Lazio Meridionale, Via Gaetano di Biasio 43, 03043 Cassino, Italy
^4 Dipartimento di Ingegneria Elettrica e delle Tecnologie dell'Informazione, Università degli Studi di Napoli Federico II, via Claudio 21, 80125 Napoli, Italy
^5 CREATE Consortium, Università degli Studi di Napoli Federico II, via Claudio 21, 80125 Napoli, Italy
^6 ICREA, Barcelona, Spain
^7 See the author list of ref. <cit.>

federico.cipolletta@bsc.es





  April 29, 2024






JOREK is an advanced non–linear simulation code for studying MHD instabilities in magnetically confined fusion plasmas and their control and/or mitigation. A free–boundary and resistive wall extension was introduced via coupling to the STARWALL and CARIDDI codes, both able to provide dense response matrices describing the electromagnetic interactions between plasma and conducting structures. For detailed CAD representations of the conducting structures and high resolutions for the plasma region, memory and computing time limitations restrict the possibility of simulating the ITER tokamak. In the present work, the Singular Value Decomposition provided by routines from the ScaLAPACK library has been successfully applied to compress some of the dense response matrices and thus optimize memory usage. This is demonstrated for simulations of Tearing Mode and Vertical Displacement Event instabilities. An outlook to future applications on large production cases and further extensions of the method are discussed.


Keywords: Magnetic Plasma Confinement, MHD, linear algebra, HPC, Numerical Methods



[
    F. Cipolletta^1, N. Schwarz^2, M. Hoelzl^2, S. Ventre^3, N. Isernia^4, G. Rubinacci^5, A. Soba^1, M. J. Mantsinen^1,6, and the JOREK Team^7
    October 2023
===============================================================================================================================================





§ INTRODUCTION



The ability to perform realistic predictive simulations of large-scale plasma instabilities in future magnetic confinement devices like ITER or DEMO is essential for their safe and efficient operation. Towards this effort, it is necessary to develop powerful multi-physics models, validate them against experiments in existing fusion devices, and optimize the models algorithmically and numerically for predictive simulations with particularly high resolutions. Only after reaching the capability to produce robust simulations of large device sizes subject to challenging plasma conditions, actual predictions can be made with sufficiently high fidelity. The present article describes efforts aiming to enable simulations that were previously not feasible, by implementing and verifying matrix compression techniques in the coupling of the JOREK non-linear MHD code <cit.> with the resistive wall codes STARWALL <cit.> and CARIDDI <cit.>. 

Accounting for resistive wall effects is indispensable for capturing magnetic field perturbations across the whole domain with full realism. Most notably, certain dangerous classes of instabilities like vertical displacement events (VDEs) <cit.> or resistive wall modes (RWMs) <cit.> would not be captured at all without these contributions. For recent applications of the JOREK code including instabilities at the plasma boundary and their control, major disruptions and their mitigation by shattered pellet injection, runaway electrons, VDEs, divertor and scrape-off layer dynamics, turbulence, and more, the reader is referred to the review article Ref. <cit.> and references therein. In the rest of this section, we briefly describe the JOREK code, its coupling to the resistive wall codes, and the aim of the work described in the rest of this article.

JOREK is a non-linear code for global simulations of instabilities in the realistic geometry of magnetic confinement fusion devices. The available physics models include reduced and full MHD fluid descriptions of the plasma with various extensions for two-fluid effects, neoclassical physics, kinetic treatment for specific particle species, etc., a variety of hybrid kinetic-fluid models, and (gyro)kinetic turbulence models. The spatial discretization is based on Bezier finite elements of arbitrary order combined with a toroidal Fourier expansion. The use of fully implicit time evolution methods avoids overly restrictive time–stepping conditions. The code uses a hybrid MPI-OpenMP parallelization.

Running JOREK standalone implies using an ideal wall boundary condition, considering all magnetic perturbations vanish at the boundary. Nevertheless, finite perturbations across the plasma boundary and the interaction of the plasma dynamics with conducting structures can be highly relevant for capturing the processes of interest. To achieve this, JOREK can leverage the calculations of external codes, STARWALL or CARIDDI <cit.>, which describe non-axisymmetric conducting structures and their interaction with the plasma. Those structures, including the walls and the coils, are discretized using a thin wall approximation in STARWALL and a volumetric description in CARIDDI, see Figure <ref> for examples. JOREK does not discretize the region outside its domain but, instead, uses a Green's functions approach. Therefore, large dense matrices are calculated in the resistive wall codes, which allow JOREK to evolve wall and coil currents in time while calculating the magnetic field at the boundary of the computational domain self-consistently. Inside JOREK, those matrices and the boundary condition calculation are distributed across MPI tasks over the available CPUs <cit.>. However, simulations with highly resolved wall models or high resolutions inside the plasma domain can be limited by memory consumption and computational costs. The goal of the developments reported in this article is to address those limitations.

 

The rest of the manuscript is organized as follows. In section <ref>, details of the implementation for the matrix compression will be given, describing the chosen techniques, the terminology adopted for expressing compression and memory savings, and the present status of the work. Then, in section <ref>, the results of the simulations validate the correct implementation of the methods and prove the first demonstrations for the application of the matrix compression techniques to two different relevant kinds of plasma instabilities. Finally, in section <ref>, the implementation and verification are summarized and an outlook to future work is given.



§ IMPLEMENTATION





 §.§ Matrix Compression



In the Linear Algebra literature, there are methods for obtaining matrix compression and controlling the accuracy of the result, with different assumptions on the structure and the features of the matrix itself. One of the less restrictive methods is based on the Singular Value Decomposition (SVD) (see e.g. ref. <cit.>). This factorization technique allows one to write a given matrix, A, with m rows and n columns, as


    A = U ·Σ· V^T,


being U and V^T the left and right orthogonal matrices, respectively, and Σ the diagonal matrix of the singular values. Although the representation of (<ref>) is not unique, one always exists with the diagonal elements of Σ in decreasing order. Denoting by k_A the rank of Σ, i.e. the number of non–zero singular values, it is then clear that U would have m rows and k_A columns, while V^T would have k_A rows and n columns, and that


    k_A≤min(m,n).


Then, instead of storing the matrix A in memory, one could store the two matrices given by (U ·Σ) and V^T, corresponding to a total of (n+m) · k_A numbers instead of m · n. Selecting one value k < k_A, thus neglecting the (k_A - k) smallest singular values, leads to compressing the required memory for storing the matrix A when


    (m+n) k < m n.


In other words, indicating with r_r the chosen rate of retained singular values, such that


    k = r_r min(m,n),


using (<ref>) in (<ref>), memory is saved provided


    r_r < max(m,n)/m+n.


The compressed SVD representation of matrix A with a rank k < k_A is often called the truncated SVD.

Additionally, it is useful to define a rate of memory gain, r_mg, of a given compressed matrix to measure the efficiency of compression in applications. Computing the percentual relative difference between the dimension of the original matrix and the compressed one, r_mg reads


    r_mg = ( 1.0 -  (m+n) r_r/max(m,n)) * 100 %.


Of course, r_mg can be negative, which is expected when the condition (<ref>) is not verified. For example, when r_r = 1.0, the matrix is not compressed but only factorized, therefore r_mg < 0. 

In terms of accuracy, finally, the Eckart–Young theorem proves that the truncated SVD is an optimal approximation of a given matrix at fixed rank when computing the error via the Frobenius norm (see <cit.> for a review and <cit.> for the original reference).



 §.§ Matrices chosen for the implementation


The coupling of JOREK to STARWALL or CARIDDI introduced in section <ref> adopts the virtual casing principle to express the interaction between the plasma and eddy currents in the wall by a set of static matrices. The reader is referred to Appendix C of <cit.> for a general description of all the geometrical matrices involved in the calculations.

In the present study, the matrices M^ey and M^ye of ref. <cit.> were chosen, which describe the mutual interactions between the wall and the plasma contained in it. Without entering into formal definitions, it is worth recalling the physical meaning of those matrices. On the one hand, the matrix M^ye relates the magnetic flux variations at the JOREK boundary to current variations in the external conductors. On the other, the matrix M^ey provides information on the tangent magnetic field to the JOREK boundary produced by the external currents when a superconducting shell at the JOREK boundary is considered.

Denoting by n_w the degrees of freedom of the wall and by n_bd the degrees of freedom at the boundary of the plasma region, M^ey is of dimension (n_bd× n_w), while M^ye is of dimension (n_w × n_bd) (see ref. <cit.> for the details).
As can be seen, both matrices have n_w · n_bd entries. In many cases, these matrices dominate the memory consumption when very detailed wall models are used in the simulations. Following the geometrical finite elements representation adopted in JOREK and recalled in section <ref>, one can write


    n_bd = n_B n_F = n_B ( 1 + 2 k_F),


being n_B the degrees of freedom of the poloidal Bezier elements, n_F the toroidal Fourier harmonics, where sine and cosine components are counted separately (as in the JOREK code), and k_F the Fourier mode number.

Therefore, one immediate way of changing the dimension n_bd of the chosen matrices, would be to increase the number of toroidal harmonics to be considered, as reported in section <ref>, even without varying n_w. This gives a fast way to test the implementation with different sizes of the matrices.



 §.§ Code Development and Optimization


The SVD is performed with the Scalable Linear Algebra PACKage (ScaLAPACK, see ref.s <cit.>), through the subroutine . As the matrices are static, the compression is only required once before the start of the JOREK simulation. The compression is performed by the newly developed  program. In detail, the  program should be executed after the production of response matrices from STARWALL or CARIDDI and before the simulation to be performed with JOREK. A simple sketch of this implementation is reported in algorithm <ref>. Here, after the singular values are determined by , the smallest values are eliminated based on the r_r of the user, and the truncated SVD is thus obtained.



For debugging purposes, it is also possible to print out re–combined versions of the factorized matrices, to print out the SVD analysis details, and to compress only the non–axisymmetric part of the matrix.

The isolation of the compression in a separate program allowed limiting the number of modifications to the original JOREK code. Nevertheless, the JOREK code had to be adapted to work with the chosen matrices in a factorized format. In references <cit.>, it can be seen that both M^ey and M^ye of section <ref> are used in algebraic computations involving matrix–vector products. Those operations were adopted to calculate the derived response matrices temporarily stored in the memory. Eliminating such derived response matrices and implementing a new generalized matrix–vector product instead save additional Random Access Memory during the code execution.

It is worth trying and estimating how much the computational cost is modified when considering factorized matrices. In particular, the matrices M^ey and M^ye are used to update the wall currents and to compute the vacuum boundary integral for obtaining the magnetic field terms (see Section 3.2 of <cit.>). There are 5 to 7 matrix–vector products per time iteration, depending on whether the previously used time–step is kept or changed. Therefore, estimating the computational cost of such a single calculation can give a general idea of the overall number of operations. Here, in general, one finds that the number of scalar multiplications is


    cost_aggr = m_A · n_A,


when the matrix A is aggregated. However, when considering the factorization of (<ref>) with a truncated SVD retaining k singular values as in (<ref>), then the overall cost is


    cost_compr = ( m_A · k ) + ( k · n_A ).


The direct comparison of the two costs provides the condition on the number of retained singular values to have computational cost savings:


    k < m_A · n_A/( m_A + n_A ).


Given that those operations happen at each time iteration, it was crucial to be sure that the factorization of the matrix does not deteriorate the overall performance. To this end, hybrid / parallelization was exploited, factorized matrices were distributed across MPI ranks, and the nested loops involving them have been adapted. Those optimizations allowed retaining the same computational time when adopting factorized but not compressed matrices or their original aggregated version.

To conclude the present section, an insight into the practical meaning of equations (<ref>), (<ref>), and (<ref>) is provided by the spectral analysis reported in figure <ref>. Given that the normalized singular values of matrix M_ye are slightly smaller than the ones of M_ey, in the adopted discretized geometrical model, the former  is expected to show better compressibility than the latter.





§ VERIFICATION AND TESTING



The JOREK code can simulate various kinds of plasma instability that may occur in a magnetically confined plasma <cit.>. To test the implementation presented in section <ref>, two such scenarios were selected, namely the Tearing Mode (TM) instability and the Vertical Displacement Event (VDE). The simulation setups and results of compression tests are shown in the following sections. Note that we used comparably small test cases here. In the future, computationally more expensive production simulations will be done after the modifications have been merged into the main code version.



 §.§ Tearing Mode Instability



The Tearing Mode (TM) instability is a well–known phenomenon in the field of MHD (see for example <cit.> for a classical review on this topic). Despite its theoretical conception being initially driven by astrophysical argumentations, the relative scenario was recovered and deeply studied in the experimental tokamak plasma physics field. This instability modifies the magnetic topology by reconnecting magnetic field lines such that the so–called magnetic islands are formed. Under certain circumstances, TMs can grow to significantly large amplitude and cause major disruptions, during which the confinement of the plasma is lost in an abrupt way causing large heat and electromagnetic loads for the device.




The test case considered here is linearly unstable to a tearing mode with dominant toroidal mode number 1 and the vacuum response matrices are calculated by the CARIDDI code. The geometrical structure adopted is similar to what is shown in figure <ref>. The axisymmetric component of the magnetic field (with mode number 0) is not included in the response matrices, such that an ideal wall boundary condition is applied for the axisymmetric component. This avoids dealing with vertical plasma instabilities, which are addressed separately in the following section.  The variation of the linear growth rate of the TM as a function of the wall resistivity is reported in figure <ref>. The nominal wall resistivity used in further studies below has a value of 0.72 × 10^-6 [ Ω· m ]. 


 
Before looking into the compression efficiency, the correctness of the implementation is tested. For this purpose, a case with n_w = 20322, n_B = 160, n_F = 5 is simulated both with the uncompressed matrices and with the factorized matrices without compression (r_r=1). The comparison of the magnetic and kinetic energies in the different toroidal harmonics shown in Figure <ref> shows perfect agreement.

The efficiency of the compression is now tested using different retention rates in the factorized matrices. In table <ref>, the values adopted for r_r are reported together with the size s of the corresponding truncated SVD matrices and the relative memory usage reduction r_mg given by (<ref>), computed based on the original size. It is worth noting that the first line corresponds to the same simulation performed with aggregated and uncompressed matrices and should give matching results with r_r = 1.0 for validation.
In the tests performed, accurate results could be obtained with retention rates down to r_r = 0.15 as shown in Figure <ref>. With smaller values of r_r the results were instead inaccurate. 







As a general statement, it is worth noting that the efficiency of the compression of a matrix is expected to grow with increasing dimensions of the matrix itself. Therefore, a series of tests with varying resolutions of the boundaries on the wall (n_w) and the poloidal section of the plasma (n_B) has been performed, keeping fixed n_F=5. These scans are performed with the STARWALL instead of the CARIDDI code for simplicity, thus adopting a geometrical structure like what is shown in figure <ref>, starting from a base resolution of n_w=20449, n_bd=400 and increasing either the resolution of the wall or the plasma independently. The results shown in figure <ref> indeed show better compressibility for larger problem sizes. However, it is also worth noting that the maximum compressibility obtainable via truncated SVD is limited by the minimum between the number of rows and columns of the matrix, as understandable from equations (<ref>) and (<ref>). In the studied cases, this limit is always given by the number of boundary elements on the plasma side.







 §.§ Vertical Displacement Event (VDE)



VDE refers to a loss of control of the vertical positioning of the magnetically confined plasma inside a tokamak (see <cit.> and reference therein for a recent report of simulations of this instability with the ITER geometry). The following loss of thermal and magnetic energy can reduce the machine's lifetime, therefore it is crucial to understand its dynamics with the help of numerical simulations. 

In references <cit.>, 

it is possible to see how such a phenomenon can be effectively set up with the STARWALL code and the AUG geometry, to be simulated with the JOREK code, taking into account initial perturbation on the coil current. Moreover, reference <cit.> reports VDE simulations to benchmark comparison between the codes M3D-C^1, NIMROD, and JOREK, with a simplified geometry. In addition, reference <cit.> contains results of simulations of VDE where the response was computed adopting CARIDDI. 
The VDE simulations of the present section closely follow the setup adopted in already published work, described in <cit.>, <cit.>. It consists of the following steps:


    
  1. a 2D axisymmetric phase, taking into account only n_F=0 harmonic until the magnetic axis reaches the vertical position Z = -0.1 m;
    
  2. a full 3D restart taking into account all the chosen n_F harmonics;    
    
  3. restart at the Thermal Quench (TQ) phase with constant perpendicular diffusivities, no particle source, and decreasing the
    time–step;
    
    those settings simplified the simulations, allowing their continuation, and alleviating overly restrictive limitations on the time–steps during TQ.





The VDE tests of the present paper were performed by adopting the AUG geometry inside CARIDDI (refer to the geometrical structure of figure <ref>) with n_w = 20322, n_B = 180, and different values of n_F (namely 1, 9, and 19), to evaluate the effect of the compression of the chosen matrices on the accuracy of the toroidal Fourier harmonics representation. Moreover, several values of r_r were initially taken into account, for both the M^ey and M^ye introduced in section <ref>. The details of the memory sizes involved in those VDE simulations are reported in table <ref>. As expected, going from the top row to the bottom one, it is possible to see that compressing the matrices at higher rates allows for saving more memory.

In all the simulated cases, the matrices factorized with r_r = 1.0 reproduced exactly the results of the original version of the matrix, thus verifying the implementation, see figure <ref>, where the time evolution of the vertical axis position is compared between the original and factorized but uncompressed matrices. 



On the other hand, all the simulations adopting r_r < 1.0 failed to produce accurate results, when compressing both the M^ey and M^ye response matrices. Although that preliminary outcome might indicate incompressibility of the matrices in the VDE test case, also the single M^ye compression was attempted with r_r = 0.4, 0.6, and 0.8 while leaving M^ey factorized with r_r = 1.0. 






Figures <ref> and <ref> show the time traces of the vertical position of the magnetic axis, the thermal energy content of the plasma, and the magnetic energies for selected harmonics when adopting n_F=19. The curves (colored in the online version) are obtained from uncompressed (indicated by standard) or compressed M^ye. The time considered in such plots spans from slightly before starting the 3D phase and up and beyond the TQ phase. In particular, figure <ref> shows that the results for the Z–coordinate of the magnetic axis and the total thermal energy obtained without compression are almost exactly matched by the ones obtained compressing the matrix M^ye, up to the rapid drop in the total thermal energy, happening during the TQ phase. Indeed, during the highly non-linear TQ, the Z position of the magnetic axis shown in figure <ref> becomes less regular. In this phase, the compressed simulations cannot accurately describe the evolution anymore. Figure <ref> confirms a similar observation in terms of the thermal energy content of the plasma: cases with compressed M^ye reproduce the dynamics of the original case up to the TQ phase, during which the accuracy has deteriorated. For the same cases, Figure <ref> compares the evolution of magnetic energies of selected toroidal components. Also here, the initially excellent agreement deteriorates during the TQ. It is worth noting, that the agreement in the higher harmonics is better than in the low harmonics.

The errors occurring during the TQ when adopting compression are likely due to the development of strong dynamical instabilities, usually able to break the magnetic confinement surfaces (see e.g. <cit.>). Furthermore, through numerical experiments, we found that the plasma evolution appears particularly sensitive to any possible asymmetry in the conducting structures around the walls, during the TQ. Addressing such loss of accuracy during the evolution at advanced times, when strong dynamics are active requires careful treatment of the compression method from the point of view of the modeled geometry. Of course, such a task is of interest for future developments and deserves detailed additional explorations.

Moreover, understanding the reason why, compressing only M^ye in VDE simulations produced accurate results while compressing M^ey, as well, produced inaccurate results could be related to the particular nature of the instability itself but requires further investigation and goes beyond the scope of the present work.



§ CONCLUSION AND OUTLOOK



In this article we explain, verify, and demonstrate a technique for compressing the dense “vacuum response matrices” provided by STARWALL or CARIDDI, which allow JOREK to self-consistently take into account the interaction of the plasma with conducting structures. The two matrices M^ey and M^ye that describe the mutual interaction between plasma and conducting structures are addressed, which typically dominate the memory consumption. The tests focus on two representative cases with moderate resolutions. 

In the simulation of a slowly growing TM instability, the axisymmetric component is kept fixed such that the response matrices only affect the non-axisymmetric components of the plasma evolution. Here, a strong reduction of the memory consumption can be obtained and an increasing compressibility is seen for both higher wall or plasma resolutions. 

The simulation of a more challenging 3D VDE behaves differently. Only the M^ye matrix could be compressed, while M^ey compression was unsuccessful. Furthermore, even at moderate compression levels, the initially excellent agreement deteriorates during the highly non-linear thermal quench phase. The reason for this can likely be attributed to the strong coupling between the different toroidal harmonics producing a very significant burst of the mode activity at the TQ phase.

The overall test results are considered promising and will be the basis for future large–scale production applications that go beyond the scope of this article.
 
Presently observed limitations provide insights regarding the future directions of this work. At first, a plain SVD application limits the compression by the smallest matrix dimension, as understandable from equation (<ref>) and visible from figure <ref>. The largest dimension of the matrices is clearly given by the degrees of freedom in the considered device's wall and grows significantly when increasing the resolution. Adopting more complex methods (see <cit.> and references therein for some examples) might allow higher compression even when adopting more detailed wall geometries. The limitations in the case of strong dynamical instabilities suggest that geometrical aspects should be taken into account in future work. Finally, the JOREK-STARWALL and JOREK-CARIDDI couplings are now restricted to reduced MHD and do not account for halo currents flowing directly between conducting structures and the plasma domain. A full MHD treatment presently under development will give rise to new response matrices that can be compressed with the techniques described in the present paper and their successors.

Despite the limitations mentioned above, the novel implementation of the response matrices compression technique via SVD presented in this manuscript provided useful results that are already applicable while paving the way toward further improvements.



§ ACKNOWLEDGMENTS


This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 — EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. Some of this work was carried out on the Marconi-Fusion supercomputer operated by CINECA. Code development and testing of the presented work have also been performed on the MareNostrum 4 cluster at the BSC, through the allocation of the Computer Applications in Science and Engineering (CASE) Department.



§ DATA AVAILABILITY


The data that support the findings of this study are available from the corresponding author upon reasonable request.





