Parallelizing TORIC

The computational problem in TORIC is one of matrix inversion. The discretization of the system of equations in Eqs. (1) leads to a block tri-diagonal system with $3\times N_r$ blocks. A weak variational formulation of the equations is used with cubic Hermite polynomials in the radial dimension, which for the FLR system, creates the sparse tri-diagonal structure. The Fourier decomposition decouples the toroidal dimension for axisymmetric devices so that it only enters parametrically as the toroidal mode number. Multiple simulations at different toroidal modes can be used to build up a complete three dimensional spectrum. The poloidal modes are coupled by the free space operators and the plasma response and create dense blocks which contain $(2 \times 3 \times N_m)^2$ elements (a factor of two for the real and imaginary parts of the three vector components.) Thus, the computational resources to solve for the fields quickly exceed the available memory of a single processor. For example, given an available memory of 2 Gigabytes of RAM, we would be limited to a maximum of approximately 150 radial elements by 128 poloidal modes, which is insufficient for the problems discussed above.

To take advantage of scalable architectures, the code has been parallelized in the power reconstruction and the matrix inversion by using the ScaLAPACK [9] library of parallelized linear algebra routines as well as direct use of the message passing interface (MPI). In this way, the limiting memory requirements of the large blocks of the block diagonal system are distributed across multiple processors, and so the problem size is limited in principal only by the available number of processors. The power and current deposition calculations are independent of flux surface, so by distributing those calculations along that dimension, those calculations are reduced from 50% to less than 2% of a typical calculation time.