Optimization of the radiation transport calculation for quasi-one-dimensional model of the ionizing gas flows

Algorithms were developed to solve the problem of radiation transport based on the method of long characteristics, taking into account its optimization in the framework of the quasi-one-dimensional approximation used to study the axisymmetric flows of the ionizing gas in coaxial channels of plasma accelerators. The optimization of calculations of the integral radiation parameters allowed to significantly reduce the calculation time in comparison with three-dimensional model of radiation transport. Calculations of the radiation energy density and the radiation energy flux density, based on various methods, have shown that new algorithms provide high-quality solution of the radiation transport problem.


Introduction
The study of radiation transport is one of the actual problems of modern computational plasma dynamics. The radiation leads to the redistribution of energy in the plasma volume, and also affects the ionization process in the channel of the quasi-stationary plasma accelerator (QSPA) [1][2][3][4][5][6][7][8][9], schematically shown in Fig. 1a. The plasma in the presence of an azimuthal magnetic field and a predominantly radial plasma current is accelerated in the QSPA due to the Ampere force c / H j , where j is the current in the plasma. Different models of radiation transport in the channel of the QSPA have been developed earlier, including the model of the radiant heat conductivity and the diffusion approach (see e.g. [10]).
Recently, the 3D model of radiation transport in the QSPA has been developed based on the method of long characteristics (see e.g. [11][12][13]), which allows to take into account the details of the geometry of the accelerator channel and the formation of shadow regions [14]. The solution of the spectral 2 radiation transport equation in three-dimensional formulation is time-consuming task in terms of computer time. In this connection, it became necessary to construct the simplified model that would allow qualitatively and quickly to find the solution of the radiation transport problem based to study the ionizing gas flow in the quasi-one-dimensional approximation, taking into account the channel geometry, the emission spectrum and the possibility of parallel implementation of computational code.
The quasi-one-dimensional approximation is used to numerically simulate the axisymmetric flows in the narrow channel [15][16][17] and allows to effectively reduce the calculation time. The study of the ionizing gas flow with the formation of an ionization front in the QSPA channel is carried out on the basis of the MHD models of various levels of complexity, using the local thermodynamic equilibrium (LTE) approximation (see e.g. [14,15]) and the modified diffusion approximation (MDA) for the equation of the ionization and recombination kinetics [16][17]. The most complete MHD model includes the set of equations of the level kinetics for populations of atomic levels. In this paper, we propose the appropriate optimization method to calculate the radiation transport in quasi-onedimensional model of observable axisymmetric flows. Calculations based on the developed algorithm are compared with the results obtained by means of the 3D model of radiation transport. The proposed method of the simplified calculation was used to study the effect of radiation transport on the ionization process in the QSPA channel.

Radiation transport equation
Here, the radiation intensity   Ω r, in the direction of the spatial angle Ω is a function of coordinate r and photon frequency  . The corresponding integrals of the radiation intensity determine the spectral density of the radiation energy () U  r and the spectral flux of the radiation energy ()  Wr: The equations of radiation magnetic gas dynamics contain the integral values of the radiation Further, we consider one-group approximation and will not indicate the index  in (1) and (2). Figure 1(b) shows the geometry of plasma accelerator channel and the grid in quasi-one-dimensional approximation. The line marked in red in the figure corresponds to the middle radius in the channel.
The integral radiation parameters will calculate at points located on the middle line between the electrodes. If parameter z N is the number of grid nodes along the z axis, equal to  (1) depend on the state of the medium, its density and temperature, as well as the photon energy h . These values are determined by means of the known relations (see e.g. [14,[18][19][20][21][22]).
Graphs of typical dependences   and   on photon energy are presented in figure 2 for hydrogen plasma taking into account the first 20 energy levels at characteristic temperature eV T 1  and particle density

Calculation of radiation field based on the method of long characteristics
The radiation transport equation (1) is the first order equation for which any ray or direction of the radiation propagation can be considered as a characteristic. The phenomenon of refraction of rays is not taken into account.
We consider a solution on the characteristics (see e.g. [11][12][13][14]). Figure 3 shows one-dimensional spatial grid, where points correspond to the nodes at which the intensity is determined in order to then find the integral radiation parameters. Within each cell, the absorption coefficient and emissivity are considered constant. Solution on the characteristic with a direction i Ω at node n is constructed from the boundary of the computational domain at b r r  with the known intensity value bi ( , ) I r Ω :  using the relations (2). In this case, the integrals correspond to the sums, in which the frequency index is not specified for convenience: The results of solving the problem of radiation transport in the 3D formulation will be presented for the points of the middle coordinate line between the electrodes, shown in figure 1(b). The calculation of the ionizing gas flow with the formation of an ionization front was carried out on the basis of the MHD model taking into account the ionization and recombination kinetics in the framework of the MDA approximation [16,17,22]. The distributions of density and temperature in the channel are presented in figure 4 for the following characteristic values of the gas parameters at the channel inlet:   In order to solve the radiation transport equation, a calculation of the large number of exponents is required, which is approximately equal to twice the number of terms in formula (3) and is proportional to the number of layers or segments of the computational domain (see figure 3). On average, each characteristic passes through 50 s N  layers. Then the total number of calculated exponents for one spectral group will be approximately equal . This explains the long time to calculate the radiation field. Ways to optimize the calculation of radiation transport are associated with a decrease in the number of the calculated exponents.

Optimization of radiation transport calculation
The number of operations in the method of long characteristics and, accordingly, the calculation time depend on the number of rays  N of the angular grid. The time needed to solve the problem of radiation transport can be significantly reduced for the quasi-one-dimensional distributions of temperature and density, taking into account the layered structure of the distribution of thermodynamic parameters in the channel.

Optimization by mean of the rays coming from boundary
Summation over the directions of the angular grid can be replaced by summation over the elements of the boundary. In this case, the sum of the exponents in formulas (3) and (4)   N is the number of ring elements of the boundary. We consider in more detail the sums over the group of rays O nb in these relations and introduce the notation:   6 We take into account that only the flow in the direction of the z axis is of interest in the quasi-onedimensional formulation of the problem. Sums (5) can be represented in more detail: Here, the summation is carried out over the cells through which the characteristic passes from the node with the index n to the cell with the number m . As a result, the sum is allocated that does not depend on the direction of the ray and cos i  .
Substituting formulas (7) into relations (6), we obtain the following dependences:  . The preparatory procedure based on the 3D model can be performed in sufficient detail for the given channel geometry. As a result, data tables are formed to calculate the function parameters with different indices i and k .

Optimization by mean of the rays coming from plasma layers
The technique presented in the previous section for the boundary of the layer can significantly reduce the calculation time. We apply it to the whole layer, introducing the designation O nl as a group of rays coming from the layer l and coming to the node n of the computational grid in accordance with figure 7(b). For the radiation energy density and the radiation energy flux density, we obtain the following relations by analogy with (10) where l N is the number of layers. We assume that the number of layers or volume elements corresponds to the number of nodes in the coordinate grid . This number is 2 times less than the number of boundary ring , taking into account the surfaces of the inner and outer electrodes of plasma accelerator. Accordingly, the calculation time for the radiation transport problem is reduced by about two times compared with optimization based on the rays coming from the boundary.
It should be noted that the optimization by means of the boundary elements involves the setting of the boundary conditions on the surface of the coaxial electrodes. At the same time, the optimization method based on integration over the layers or the corresponding volume elements using relations (11) excludes the specification of boundary conditions on the surface of accelerator electrodes and involves the specification of boundary conditions at the channel inlet and outlet, which also provides the qualitative solution of the radiation transport problem.

Specification of function parameters to optimize calculations
In this section, we consider the function-parameters p ( ) ik  and r ( ) ik  , which include relations (10) and (11). The characteristic dependences of the function parameters on the optical length  are presented in figure 8.

Results of radiation transport calculations based on optimization methods
Calculations of the radiation energy density and the radiation energy flux density in the quasi-onedimensional model of the ionizing gas flow are based on the methods of summation over the elements of the boundary and summation over layers or elements of the volume. Comparison of these calculations among themselves showed that both of these optimization methods give almost identical results. The calculation time of the integral characteristics of radiation based on the method of summation over layers was reduced to 4.5 seconds. This is two times less than the calculation using the summing over the elements of the boundary, for which the calculation time was 9 seconds. In both cases, the computation time with optimization is significantly less than the computation time of the integral radiation parameters using the characteristics method in the 3D model. Results of calculations based on the standard method of characteristics were represented by solid curves in Fig. 5. In this case, the calculation time was 670 seconds for the angular grid with the number of rays 3280   N . The distributions of the radiation energy density and the radiation energy flux density calculated on the basis of the summing method over layers are represented by the dashed curves in figures 5(a) and 5(b). It can be seen that calculations based on the method of summation over layers and the standard method of characteristics give practically identical results. Thus, the use of optimization methods can significantly reduce the computation time, ensuring the quality and accuracy of solving the problem of radiation transport.
The study of radiation transport in streams of ionizing gas showed that radiation, coming from the ionization front, in certain frequency range can penetrate deep into the region of the incoming neutral gas. Therefore, along with the redistribution of the internal energy of the medium and temperature, the radiation can lead to the effect of pre-ionization in the gas flow ahead of the ionization front. The frequency in the radiation transport equation (1) plays the role of an independent parameter, which allowed to develop parallel program codes that were created using of the OpenMP technology.

Conclusion
Numerical algorithms have been developed for solving the problem of radiation transport in the axisymmetric ionizing gas flows with the formation of an ionization front in the coaxial channels of plasma accelerators. The proposed optimization methods to calculate of radiation transport in the framework of the quasi-one-dimensional MHD model of flows can significantly reduce the computation time of the integral radiation parameters in comparison with the 3D model of radiation transport based on the method of long characteristics. Comparison of the calculation results for the radiation energy density and the radiation energy flux density showed that optimization of calculations based on the summation of radiation from the ring elements of the boundary and the summation over the plasma layers allows to obtain the high-quality solutions of the problem.