Parallelization of Non-Equilibrium Radiation Transport Code

Igor E. Golovkin^{∗}  Roberto C. Mancini^{∗}  Frederick C. Harris, Jr.^{†}

Radiation transport effects on X-ray line emission are important in modeling spectroscopic-quality, synthetic spectra of mid- and high-Z tracer elements in hot dense plasmas. Plasmas produced in the laboratory are usually in non-equilibrium, i.e. the ionization balance and distribution of atomic-level populations are determined by a set of collisional and radiative atomic processes. In these conditions radiation in the plasma and the level-population distribution are interdependent and have to be self-consistently determined. This involves the simultaneous solution of a set of atomic kinetics rate equations and the radiation transport equation, a problem which is non-linear and non-local. This results in an integro-differential problem that in general can not be solved analytically [1].

The problem of non-equilibrium radiation transport is quite computationally intensive. Within the last 15 years Lambda operator techniques have been introduced and have produced robust iterative schemes [2]. We have utilized a combination of linearization and the Lambda operator approach. Our model focuses on a spectral range that covers line transitions relevant for spectroscopy diagnostics. We assume a plane parallel slab geometry with no incident radiation on either side of the slab.

The basic outline of this iterative method can be described as follows: first we introduce spatial discretization. Then for each spatial zone we: (1) linearize the atomic kinetic rate equations; (2) replace one of the equations by a particle conservation condition, since the set of atomic kinetic rate equations is redundant; (3) express radiation dependent terms in these equations through population numbers using Lambda operator; this is an integral operator based on a formal solution of the radiation transport equation; (4) solve the resulting system for corrections to population numbers; (5) update the values of the population numbers and the radiation field and then repeat all these steps until a convergence criterion is satisfied.

This procedure allows us to obtain the solution, i.e. the population numbers and radiation field for each spatial zone. This method is robust and stable, does not require a lot of memory, and it converges rapidly. The only disadvantage is that it is very computationally intensive, mainly because of the complexity of the Lambda operator.

Several approximations and acceleration techniques have been developed to reduce the computational time required for calculating the Lambda operator. Many of these methods have been shown to reduce computational time significantly; unfortunately, these methods also have several drawbacks. The first disadvantage is that the choice of an approximate Lambda operator may depend on the physical problem under consideration. This hurts the universality of the code. The second disadvantage is that, an approximate operator may contain optimization parameters that can not be known a priori but have to be found by trial and error. Finally, a situation may arise in which by reducing computational time per iteration one may have to pay the price of increasing the number of iteration steps. This may result in a small improvement, no improvement, or worse yet, a longer running time.

∗Department of Physics, University of Nevada, Reno, golovkin@physics.unr.edu; rcman@physics.unr.edu †Department of Computer Science, University of Nevada, Reno, fredh@cs.unr.edu

Because we wanted to keep our algorithm general, straightforward, and therefore easy to modify, instead of trying the accelerations we attempted to exploit the power of parallelization. This algorithm provides a very good opportunity for distributing tasks among processing elements. Each iteration step consist of two steps: building the system of equations Ax=b (i.e. setting up matrix A and vector b), and solving it. The matrix’s dimension is equal to the number of atomic levels times the number of discretization points in space. The number of levels is determined by the physics of the problem and in many cases is not very high. The number of spatial zones can also be kept small. The matrix is not singular and therefore the system can be solved quite easily. What takes most of the computational time is setting up the matrix and the right-hand-side vector. Each non-zero element in the matrix contains a combination of triple integrals arising from the Lambda operator. These are integrals in space, angle and frequency approximated by Gauss-Legendre quadrature formulas. The first two sums can have a moderate number of terms, but for the frequency integral we must use a fine grid to get high quality spectra.

From the physics point of view each row of the matrix shows how the population of each particular atomic level in each spatial zone can be effected by all the other levels. The rows of this matrix are completely independent and thus can be calculated by a separate processor. Hence our main efforts have been focused on parallelization of the matrix set-up because an almost perfect speed-up can be expected as long as the matrix dimension is divisible by the number of processors, and each processor can be assigned the same amount of work.

We started our work by writing a sequential version of the program. This was done for two reasons. First we needed a working sequential code to serve as a diagnostic tool for our parallel code, and second, we used the sequential code as a springboard for our parallel development. The parallel code was developed and tested on our SGI Power Challenge machine with 8 processors and shared memory. Close to linear speed-up has been achieved as is shown in Figure 1.

References

1 2 3 4 5 6 1 2 3 4 5 6 Speedup  Processors ◇ ◇ ◇ ◇  ◇ Fig. 1. Problem Speed-up vs. the Number of Processors Used [1] D. Mihalas. Stellar Atmospheres. Freeman, San Francisco, CS, 2nd edition, 1978.

[2] G.B. Scharmer and M. Carlsson. A new approach to multi-level non-lte radiative transfer

problems. Journal of Computational Physics, 59:56–80, 1985.