Issue 
A&A
Volume 507, Number 3, December I 2009



Page(s)  1815  1818  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/200912491  
Published online  01 October 2009 
A&A 507, 18151818 (2009)
A conjugate gradient method for solving the nonLTE line
radiation transfer problem
(Research Note)
F. Paletou  E. Anterrieu
Laboratoire d'Astrophysique de ToulouseTarbes, Université de Toulouse, CNRS, 14 Av. E. Belin, 31400 Toulouse, France
Received 14 May 2009 / Accepted 26 August 2009
Abstract
This study concerns the fast and accurate solution of the line
radiation transfer problem, under nonLTE conditions. We propose and
evaluate an alternative iterative scheme to the classical ALIJacobi
method, and to the more recently proposed GaussSeidel and successive
overrelaxation (GS/SOR) schemes. Our study is indeed based on applying
a preconditioned biconjugate gradient method (BiCGP). Standard tests,
in 1D plane parallel geometry and in the frame of the twolevel
atom model with monochromatic scattering are discussed. Rates of
convergence between the previously mentioned iterative schemes are
compared, as are their respective timing properties. The smoothing
capability of the BiCGP method is also demonstrated.
Key words: radiative transfer  methods: numerical
1 Introduction
The solution of the radiative transfer equation, under nonLTE conditions, is a classical problem in astrophysics. Many numerical methods have been used since the beginning of the computer era, from difference equations methods (e.g., Feautrier 1964; Cuny 1967; Auer & Mihalas 1969) to iterative schemes (e.g., Cannon 1973; Scharmer 1981; Olson et al. 1986).
Since the seminal paper of Olson et al. (1986), the socalled ALI or Accelerated Iteration scheme is nowadays one of the most popular methods for solving complex radiation transfer problems. When using the diagonal of the operator as approximate operator (see e.g., Mihalas 1978), it is a Jacobi iterative scheme. It has been generalized for multilevel atom problems (Rybicki & Hummer 1991), for multidimensional geometries (Auer & Paletou 1994; Auer et al. 1994; van Noort et al. 2002), and for polarized radiation transfer (e.g., Faurobert et al. 1997; Trujillo Bueno & Manso Sainz 1999).
Despite the many success of the ALI method, Trujillo Bueno & Fabiani Bendicho (1995) proposed a novel iterative scheme for the solution of the nonLTE radiation transfer. They adapted GaussSeidel and successive overrelaxation (GS/SOR) iterative schemes, wellknown in applied mathematics as superior, in terms of convergence rate, to the ALIJacobi iterative scheme. It is interesting to note that, besides their application to radiation transfer, both Jacobi and GS/SOR iterative schemes were proposed during the 19th century.
Conjugate gradientlike, hereafter CGlike, iterative methods were proposed later by Hestenes & Stiefel (1952). They are distinct than the ALIJacobi and GS/SOR methods. Unlike the last two, CGlike methods are socalled nonstationary iterative methods (see e.g., Saad 2008). Very few articles discuss its application to the radiation transfer equation (see e.g., Amosov & Dmitriev 2005) and, to the best of our knowledge, they have not been properly evaluated, so far, for astrophysical purposes.
Such an evaluation is the scope of the present research note. This alternative approach is relevant to the quest for ever faster and accurate numerical methods for solving the nonLTE line radiation transfer problem.
2 The iterative scheme
In the twolevel atom case and with complete frequency
redistribution, the nonLTE line source function is usually written as
where is the optical depth, the collisional destruction probability^{}, B the Planck function and the usual mean intensity:
where the optical depth dependence of the specific intensity (and, possibly, of the line absorption profile ) is omitted for simplicity.
The mean intensity is usually written as the formal solution of the
radiative transfer equation
(3) 
or, in other terms, the solution of the radiation transfer equation for a known source function (see e.g., Mihalas 1978).
Let us define
,
the unknown
and the right hand side
.
Then, the
solution to the radiative transfer equation is equivalent to solving
the system of equations:
Ax = b .  (4) 
Among the various forms of CGbased methods, the biconjugate gradient (hereafter BiCG) method is more appropriate for the case of operators A, which are not symmetric positive definite. This is indeed the case for our radiative transfer problem, as illustrated in Fig. 1 where we display the structure of the operator A for the case of a symmetric 1D grid of maximum optical depth , discretized with 8 spatial points per decade and using .
Figure 1: Typical structure of the A operator, computed for the case of a selfemitting 1D, plane parallel finite slab. Such an operator is, more generally, nonsymmetric positive definite. For the sake of getting a better contrast for the offdiagonal parts, we magnified the negative values of the original operator by a factor of 5, and we display in fact the opposite of the modified operator. 

Open with DEXTER 
A general description of the BiCG gradient method can be found in several classical textbooks of applied mathematics (see e.g., Saad 2008; it is important to note that this method requires the use of the transpose of the A operator, ). However, the BiCG method remains efficient until the A operator becomes very badly conditioned, i.e. for very high albedo cases, typically when . For these cases, it is therefore crucial to switch to a more efficient scheme, using preconditioning of the system of equations.
Hereafter, we derive the algorithm for the BiCG method with
preconditioning (hereafter BiCGP). In such a case, one instead seeks
for a solution of
M^{1} A x = M^{1} b ,  (5) 
where the matrix M should be easy to invert. A ``natural'' choice for the nonLTE radiation transfer problem is to precondition the system of equations with the diagonal of the full operator A.
From an initial guess for the source function,
x_{0}=S^{(0)}, we compute a
residual r_{0} such as
r_{0} = b  Ax_{0} .  (6) 
We then set the second residual such that , where (a,b) means the inner product between vectors a and b. For all the cases considered by us, setting proved to be adequate.
Now, the BiCGP algorithm consists in running the following iterative
scheme until convergence, where j is the iterative scheme index. The
first steps consist in solving
and
(8) 
These steps are both simple and fast to compute, since M is a diagonal operator. Then one computes
(9) 
If , then the method fails; otherwise, the algorithm continues with the following computations. If j=1, then set p_{j}=z_{j1} and . For j>1, compute
(10) 
(11) 
and
(12) 
Then make q_{j}=Ap_{j} and , as well as,
(13) 
Finally, one advances the source function according to
(14) 
and the two residuals:
(15) 
and,
This process, from Eqs. (7) to (16), is then repeated to convergence. We find that a convenient stopping criterion is to interrupt the iterative process when
Indeed, this stopping criterion guarantees that the floor value of the true error, defined in Sect. 3.1, is always reached.
Figure 2: Solutions for the monochromatic scattering in an isothermal atmosphere, obtained for different values of , ranging from 10^{4} to 10^{12}. Both the expected surface law and thermalization lengths are perfectly recovered. 

Open with DEXTER 
3 Results
The solution of the nonLTE radiative transfer problem for a constant property, planeparallel 1D slab is a standard test for any new numerical method. Indeed, the socalled monochromatic scattering case allows direct comparison of the numerical solutions to the analytical ``solution'', , which can be derived in Eddington's approximation (see e.g., the discussion in Sect. 5. of Chevallier et al. 2003).
3.1 Scaling laws for nonLTE source function solutions
In Fig. 2, we display the run of the source function with optical depth, for a selfemitting slab of total optical thickness , using a 9point per decade spatial grid, a onepoint angular quadrature with , and values ranging from 10^{4} to 10^{12}. Both the surface law and thermalization lengths are perfectly recovered, even for the numerically difficult cases of large albedos ( ).
The ``true'' error, that is,
(18) 
is displayed in Fig. 3, respectively for the GS, SOR, and BiCGP iterative schemes. SOR was used with a overrelaxation parameter, and in that case.
Figure 3: History of the true error vs. iteration number for , respectively, with GaussSeidel, SOR (with ), and BiCGP with the diagonal of A as a preconditioner. 

Open with DEXTER 
3.2 Timing properties
The typical timing properties of the BiCGP scheme, depending on the optical depth grid refinement expressed as the number of point per decade, , for a , and selfemitting slab, are summarized in Table 1. While each GS/SOR iteration is about 20% longer than one ALIiteration, the increase in time per iteration, with the number of depth points, for BiCGP is significant. However, the balance between computing time and reduction in the number of iterations always favors the BiCGP scheme vs. GS/SOR and, a fortiori, ALI.
Table 1: Timing properties of the BiCGP iterative scheme.
3.3 Sensitivity to the spatial refinement
We adopted the same slab model as the one reported in the previous section, with , and we allowed the grid refinement to vary from 5 to 20 points per decade, with a step of 3.
Figure 4: Evolution of the true error for varying spatial grid refinements, from 5 to 20 points per decade, with steps of 3. The value of the true error plateau decreases with increasing grid refinement. 

Open with DEXTER 
3.4 Smoothing capability
A strong argument for the GaussSeidel iterative scheme resides in its smoothing capability. This point is particularly important in the framework of multigrid methods for the solution of nonLTE radiative transfer problems (see e.g., Fabiani Bendicho et al. 1997).
Figure 5: Initialization and evolution of the error, , during the 5 first iterations of the BiCGP scheme. The highfrequency component introduced ab initio in the source function is removed after the very first iteration of BiCGP. Other curves, of continuously decreasing amplitudes, correspond to iterations 2 to 5. 

Open with DEXTER 
After the very first iterate of the BiCGP iterative process, the error has an amplitude of 0.75 and, more important, it is already completely smoothed. Other error curves displayed in Fig. 5, of continuously decreasing amplitudes, are the successive errors at iterates 2 to 5.
This simple numerical experiment demonstrates that BiCGP can also be used as an efficient smoother for multigrid methods.
4 Conclusion
We propose an alternative method for solving the nonLTE radiation transfer problem. Preliminary but standards tests for the twolevel atom case in 1D plane parallel geometry have been successful. In such a case, the timing properties of the preconditioned BiCG method are comparable to the ones of ALI, GS, and SOR stationnary iterative methods. However, the convergence rate of BiCGP is significantly better.
The main potential drawback of the BiCGP method is in using the operator. A detailed tradeoff analysis between overcomputing time and an improved convergence rate, with respectto the ones of GS/SOR iterative schemes, is currently being conducted.
Multilevel and multidimensional cases will be considered further. However, the use of BiCGP does not require the cumbersome modifications of multidimensional formal solvers required by GS/SOR methods (see e.g., Léger et al. 2007).
AcknowledgementsWe are grateful to Drs. Bernard Rutily and Loïc Chevallier for numerous discussions about how to solve the radiative transfer equation.
References
 Amosov, A. A., & Dmitriev, V. V. 2005, MPEI Bulletin, 6, 5 (in russian)
 Auer, L. H., & Mihalas, D. 1969, ApJ, 158, 641 [NASA ADS] [CrossRef]
 Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS]
 Auer, L. H., Fabiani Bendicho, P., & Trujillo Bueno, J. 1994, A&A, 292, 599 [NASA ADS]
 Cannon, C. J. 1973, ApJ, 185, 621 [NASA ADS] [CrossRef]
 Chevallier, L., Paletou, F., & Rutily, B. 2003, A&A, 411, 221 [NASA ADS] [CrossRef] [EDP Sciences]
 Cuny, Y. 1967, Ann. Ap., 30, 143 [NASA ADS]
 Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. H. 1997, A&A, 324, 161 [NASA ADS]
 Faurobert, M., Frisch, H., & Nagendra, K. N. 1997, A&A, 322, 896 [NASA ADS]
 Feautrier, P. 1964, C.R.A.S., 258, 3189
 Hestenes, M. R., & Stiefel, E. 1952, Journal of Research of the National Bureau of Standards, 49(6), 409
 Léger, L., Chevallier, L., & Paletou, F. 2007, A&A, 470, 1 [NASA ADS] [CrossRef] [EDP Sciences]
 Mihalas, D. 1978, Stellar Atmospheres (San Francisco: Freeman)
 Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, JQSRT, 35, 431 [NASA ADS]
 Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171 [NASA ADS]
 Scharmer, G. B. 1981, ApJ, 249, 720 [NASA ADS] [CrossRef]
 Saad, Y. 2008, Iterative methods for sparse linear systems (Philadelphia: SIAM)
 Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646 [NASA ADS] [CrossRef]
 Trujillo Bueno, J., & Manso Sainz, R. 1999, ApJ, 516, 436 [NASA ADS] [CrossRef]
 van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066 [NASA ADS] [CrossRef]
Footnotes
 ... probability^{}
 The albedo is more commonly used in general studies of the radiation transfer equation.
All Tables
Table 1: Timing properties of the BiCGP iterative scheme.
All Figures
Figure 1: Typical structure of the A operator, computed for the case of a selfemitting 1D, plane parallel finite slab. Such an operator is, more generally, nonsymmetric positive definite. For the sake of getting a better contrast for the offdiagonal parts, we magnified the negative values of the original operator by a factor of 5, and we display in fact the opposite of the modified operator. 

Open with DEXTER  
In the text 
Figure 2: Solutions for the monochromatic scattering in an isothermal atmosphere, obtained for different values of , ranging from 10^{4} to 10^{12}. Both the expected surface law and thermalization lengths are perfectly recovered. 

Open with DEXTER  
In the text 
Figure 3: History of the true error vs. iteration number for , respectively, with GaussSeidel, SOR (with ), and BiCGP with the diagonal of A as a preconditioner. 

Open with DEXTER  
In the text 
Figure 4: Evolution of the true error for varying spatial grid refinements, from 5 to 20 points per decade, with steps of 3. The value of the true error plateau decreases with increasing grid refinement. 

Open with DEXTER  
In the text 
Figure 5: Initialization and evolution of the error, , during the 5 first iterations of the BiCGP scheme. The highfrequency component introduced ab initio in the source function is removed after the very first iteration of BiCGP. Other curves, of continuously decreasing amplitudes, correspond to iterations 2 to 5. 

Open with DEXTER  
In the text 
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.