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/0004-6361/200912491 | |
Published online | 01 October 2009 |
A&A 507, 1815-1818 (2009)
A conjugate gradient method for solving the non-LTE line
radiation transfer problem
(Research Note)
F. Paletou - E. Anterrieu
Laboratoire d'Astrophysique de Toulouse-Tarbes, 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 non-LTE conditions. We propose and
evaluate an alternative iterative scheme to the classical ALI-Jacobi
method, and to the more recently proposed Gauss-Seidel and successive
over-relaxation (GS/SOR) schemes. Our study is indeed based on applying
a preconditioned bi-conjugate gradient method (BiCG-P). Standard tests,
in 1D plane parallel geometry and in the frame of the two-level
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 BiCG-P method is also demonstrated.
Key words: radiative transfer - methods: numerical
1 Introduction
The solution of the radiative transfer equation, under non-LTE 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 so-called 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 multi-dimensional 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 non-LTE radiation transfer. They adapted Gauss-Seidel and successive over-relaxation (GS/SOR) iterative schemes, well-known in applied mathematics as superior, in terms of convergence rate, to the ALI-Jacobi 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 gradient-like, hereafter CG-like, iterative methods were proposed later by Hestenes & Stiefel (1952). They are distinct than the ALI-Jacobi and GS/SOR methods. Unlike the last two, CG-like methods are so-called non-stationary 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 non-LTE line radiation transfer problem.
2 The iterative scheme
In the two-level atom case and with complete frequency
redistribution, the non-LTE line source function is usually written as
where


![[*]](/icons/foot_motif.png)

where the optical depth dependence of the specific intensity


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 CG-based methods, the bi-conjugate 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


![]() |
Figure 1: Typical structure of the A operator, computed for the case of a self-emitting 1D, plane parallel finite slab. Such an operator is, more generally, non-symmetric positive definite. For the sake of getting a better contrast for the off-diagonal 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 BiCG-P). 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 non-LTE 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,
x0=S(0), we compute a
residual r0 such as
r0 = b - Ax0 . | (6) |
We then set the second residual



Now, the BiCG-P 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


![]() |
(10) |
![]() |
(11) |
and
![]() |
(12) |
Then make qj=Apj and

![]() |
(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
|
Open with DEXTER |
3 Results
The solution of the non-LTE radiative transfer problem for a constant
property, plane-parallel 1D slab is a standard test for any new
numerical method. Indeed, the so-called 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 non-LTE source function solutions
In Fig. 2, we display the run of the source function with
optical depth, for a self-emitting slab of total optical
thickness
,
using a 9-point per decade
spatial grid, a one-point 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 BiCG-P iterative schemes. SOR was used with a


![]() |
Figure 3:
History of the true error vs. iteration number for
|
Open with DEXTER |
3.2 Timing properties
The typical timing properties of the BiCG-P scheme, depending on the
optical depth grid refinement expressed as the number of point per
decade, ,
for a
,
and
self-emitting slab, are summarized in Table 1. While each GS/SOR iteration is about 20% longer than one
ALI-iteration, the increase in time per iteration, with the number of
depth points, for BiCG-P is significant. However, the balance between
computing time and reduction in the number of iterations always favors
the BiCG-P scheme vs. GS/SOR and, a fortiori, ALI.
Table 1: Timing properties of the BiCG-P 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 Gauss-Seidel iterative scheme resides in its smoothing capability. This point is particularly important in the framework of multi-grid methods for the solution of non-LTE radiative transfer problems (see e.g., Fabiani Bendicho et al. 1997).
![]() |
Figure 5:
Initialization and evolution of the error,
|
Open with DEXTER |


After the very first iterate of the BiCG-P 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 BiCG-P can also be used as an efficient smoother for multi-grid methods.
4 Conclusion
We propose an alternative method for solving the non-LTE radiation transfer problem. Preliminary but standards tests for the two-level 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 BiCG-P is significantly better.
The main potential drawback of the BiCG-P method is in using
the
operator. A detailed trade-off analysis between
over-computing time and an improved convergence rate, with respectto
the ones of GS/SOR iterative schemes, is currently being
conducted.
Multi-level and multi-dimensional cases will be considered further. However, the use of BiCG-P does not require the cumbersome modifications of multi-dimensional 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 BiCG-P iterative scheme.
All Figures
![]() |
Figure 1: Typical structure of the A operator, computed for the case of a self-emitting 1D, plane parallel finite slab. Such an operator is, more generally, non-symmetric positive definite. For the sake of getting a better contrast for the off-diagonal 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
|
Open with DEXTER | |
In the text |
![]() |
Figure 3:
History of the true error vs. iteration number for
|
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,
|
Open with DEXTER | |
In the text |
Copyright ESO 2009
Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.