A conjugate gradient method for solving the non-LTE line radiation transfer problem

F. Paletou; E. Anterrieu

doi:10.1051/0004-6361/200912491

Home

All issues

Volume 507 / No 3 (December I 2009)

A&A, 507 3 (2009) 1815-1818

Full HTML

Free Access

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 $\Lambda$ -Iteration scheme is nowadays one of the most popular methods for solving complex radiation transfer problems. When using the diagonal of the $\Lambda$ 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

$\begin{displaymath}S(\tau) = (1 - \varepsilon) \bar{J}(\tau) + \varepsilon B(\tau) , \end{displaymath}$

(1)

where $\tau$ is the optical depth, $\varepsilon$ the collisional destruction probability, B the Planck function and $\bar{J}$ the usual mean intensity:

$\begin{displaymath}\bar{J} = \oint \displaystyle{{{{\rm d} \Omega} \over {4 \pi}... ...nt_{-\infty}^{+\infty} {\phi_{\nu} I_{\nu \Omega} \rm d \nu} , \end{displaymath}$

(2)

where the optical depth dependence of the specific intensity $I_{\nu \Omega}$ (and, possibly, of the line absorption profile $\phi_{\nu}$ ) is omitted for simplicity.

The mean intensity is usually written as the formal solution of the radiative transfer equation

$\begin{displaymath}\bar{J} = \Lambda [S] , \end{displaymath}$

(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 $A \equiv {Id} - (1-\varepsilon)\Lambda$ , the unknown $x \equiv S$ and the right hand side $b \equiv \varepsilon B$ . 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 $\tau_{\rm max}=2 \times 10^{8}$ , discretized with 8 spatial points per decade and using $\varepsilon=10^{-6}$ .

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{12491fg1.eps} \end{figure}$	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, $A^{\rm T}$ ). However, the BiCG method remains efficient until the A operator becomes very badly conditioned, i.e. for very high albedo cases, typically when $\varepsilon < 10^{-6}$ . 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, x₀=S⁽⁰⁾, we compute a residual r₀ such as

r₀ = b - Ax₀ .

(6)

We then set the second residual $\tilde{r}_{0}$ such that $(r_{0},\tilde{r}_{0}) \neq 0$ , where (a,b) means the inner product between vectors a and b. For all the cases considered by us, setting $r_{0}=\tilde{r}_{0}$ proved to be adequate.

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

M z_j-1 = r_j-1 ,

(7)

and

$\begin{displaymath}M^{\rm T} \tilde{z}_{j-1} = \tilde{r}_{j-1} . \end{displaymath}$

(8)

These steps are both simple and fast to compute, since M is a diagonal operator. Then one computes

$\begin{displaymath}\rho_{j-1} = \left(z^{\rm T}_{j-1},\tilde{r}_{j-1}\right) . \end{displaymath}$

(9)

If $\rho_{j-1}=0$ , then the method fails; otherwise, the algorithm continues with the following computations. If j=1, then set p_j=z_j-1 and $\tilde{p}_{j}=\tilde{z}_{j-1}$ . For j>1, compute

$\begin{displaymath}\beta_{j-1} = \rho_{j-1}/\rho_{j-2} , \end{displaymath}$

(10)

$\begin{displaymath}p_{j} = z_{j-1} + \beta_{j-1} p_{j-1} , \end{displaymath}$

(11)

and

$\begin{displaymath}\tilde{p}_{j} = \tilde{z}_{j-1} + \beta_{j-1} \tilde{p}_{j-1} . \end{displaymath}$

(12)

Then make q_j=Ap_j and $\tilde{q}_{j}=A^{\rm T} \tilde{p}_{j}$ , as well as,

$\begin{displaymath}\alpha_{j} = { {\rho_{j-1}} \over {({\tilde{p}^{\rm T}_{j}},{q_{j}})} } \cdot \end{displaymath}$

(13)

Finally, one advances the source function according to

$\begin{displaymath}x_{j} = x_{j-1} + \alpha_{j} p_{j} , \end{displaymath}$

(14)

and the two residuals:

$\begin{displaymath}r_{j} = r_{j-1} - \alpha_{j} q_{j} , \end{displaymath}$

(15)

and,

$\begin{displaymath}\tilde{r}_{j} = \tilde{r}_{j-1} - \alpha_{j} \tilde{q}_{j} . \end{displaymath}$

(16)

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

$\begin{displaymath}{ {\Vert r \Vert _{2}} \over {\Vert b \Vert _{2}}} < 10^{-2} . \end{displaymath}$

(17)

Indeed, this stopping criterion guarantees that the floor value of the true error, defined in Sect. 3.1, is always reached.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{12491fg2.eps} \end{figure}$	Figure 2: Solutions for the monochromatic scattering in an isothermal atmosphere, obtained for different values of $\varepsilon$ , ranging from 10^-4 to 10^-12. Both the expected $\sqrt {\varepsilon }$ surface law and thermalization lengths $\approx$ $1/\sqrt {\varepsilon }$ are perfectly recovered.
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'', $S_{\rm Edd}$ , 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 $\tau_{\rm max}=2 \times 10^{8}$ , using a 9-point per decade spatial grid, a one-point angular quadrature with $\mu=\pm 1/\sqrt{3}$ , and $\varepsilon$ values ranging from 10^-4 to 10^-12. Both the $\sqrt {\varepsilon }$ surface law and thermalization lengths $\approx 1/\sqrt{\varepsilon}$ are perfectly recovered, even for the numerically difficult cases of large albedos ( $a \lessapprox 1$ ).

The ``true'' error, that is,

$\begin{displaymath}{T}_{\rm e} = {{\rm max}} \left( { {\mid S^{(j)} - S_{{\rm Edd}} \mid} \over {S_{\rm Edd}} } \right) , \end{displaymath}$

(18)

is displayed in Fig. 3, respectively for the GS, SOR, and BiCG-P iterative schemes. SOR was used with a $\omega =1.5$ over-relaxation parameter, and $\varepsilon =10^{-12}$ in that case.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{12491fg3.eps} \end{figure}$	Figure 3: History of the true error vs. iteration number for $\varepsilon =10^{-12}$ , respectively, with Gauss-Seidel, SOR (with $\omega =1.5$ ), and BiCG-P with the diagonal of A as a preconditioner.
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, $N_{\tau}$ , for a $\tau_{\rm max}=2 \times 10^{8}$ , and $\varepsilon=10^{-8}$ 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 $\varepsilon=10^{-8}$ , and we allowed the grid refinement to vary from 5 to 20 points per decade, with a step of 3.

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{12491fg4.eps} \end{figure}$	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

Figure 4 shows the sensitivity of the method with grid refinement. Increasing grid refinement corresponds to a decrease in the floor value of the true error, which demonstrates the need for a sufficient spatial resolution for the sake of accuracy of the numerical solutions. The number of iterative steps needed to fulfill the stopping criterion defined in Eq. (17) is practically linear in the number of depth points.

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).

$\begin{figure} \par\includegraphics[width=8.8cm,clip]{12491fg5.eps} \end{figure}$	Figure 5: Initialization and evolution of the error, $S^{(j)}-S^{(\infty )}$ , during the 5 first iterations of the BiCG-P scheme. The high-frequency component introduced ab initio in the source function is removed after the very first iteration of BiCG-P. Other curves, of continuously decreasing amplitudes, correspond to iterations 2 to 5.
Open with DEXTER

We have reproduced the test initially proposed by Trujillo Bueno & Fabiani Bendicho (1995, see their Fig. 9) for the GS method. Here we consider a semi-infinite slab of maximum optical thickness 10⁵, sampled with a 9-point per decade grid, and $\varepsilon=10^{-6}$ . The test consists in ab initio injecting a high-frequency component into the error, $S^{(0)}-S^{(\infty)}$ , on the initial guess for the source function. In Fig. 5, the initial error we used is easily recognized as the hyperbolic tangent-like curve, to which a high-frequency component is added.

After the very first iterate of the BiCG-P iterative process, the error has an amplitude of $\approx$ 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 $A^{\rm T}$ 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).

Acknowledgements

We 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 $a=(1-\varepsilon)$ is more commonly used in general studies of the radiation transfer equation.

All Tables