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

F. Paletou, Eric Anterrieu

► To cite this version:

F. Paletou, Eric Anterrieu. A conjugate gradient method for solving the non-LTE line radiation transfer problem. Astronomy and Astrophysics - A&A, 2009, 507 (3), pp.1815-1818. 10.1051/0004-6361/200912491. hal-02105965

HAL Id: hal-02105965

https://hal.science/hal-02105965

Submitted on 20 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A conjugate gradient method for solving the non-LTE line radiation transfer problem
*(Research Note)*

F. Paletou and E. Anterrieu

Laboratoire d'Astrophysique de Toulouse-Tarbes, Université de Toulouse, CNRS, 14 Av. E. Belin, 31400 Toulouse, France
e-mail: fpaletou@ast.obs-mip.fr

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

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

where $\tau$ is the optical depth, $\varepsilon$ the collisional destruction probability<fnref n="1" />, $B$ the Planck function and $\bar{J}$ the usual mean intensity:

$$\bar{J} = \int \frac{d\Omega}{4\pi} \int_{-\infty}^{+\infty} \phi_v I_{v\Omega} dv, \quad (2)$$

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

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

$$\bar{J} = \Lambda[S], \quad (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 \text{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. \quad (4)$$

<footnote n="1">¹ The albedo $a = (1 - \varepsilon)$ is more commonly used in general studies of the radiation transfer equation.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 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.</figcaption></figure>

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_{\max} = 2 \times 10^8$, discretized with 8 spatial points per decade and using $\epsilon = 10^{-6}$.

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^T$). However, the BiCG method remains efficient until the A operator becomes very badly conditioned, i.e. for very high albedo cases, typically when $\epsilon < 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}Ax = M^{-1}b, \qquad (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_0 = S^{(0)}$, we compute a *residual* $r_0$ such as

$$r_0 = b - Ax_0. \qquad (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

$$Mz_{j-1} = r_{j-1}, \qquad (7)$$

and

$$M^T \tilde{z}_{j-1} = \tilde{r}_{j-1}. \qquad (8)$$

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

$$\rho_{j-1} = (z_{j-1}^T, \tilde{r}_{j-1}). \qquad (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

$$\beta_{j-1} = \rho_{j-1}/\rho_{j-2}, \qquad (10)$$

$$p_j = z_{j-1} + \beta_{j-1}p_{j-1}, \qquad (11)$$

and

$$\tilde{p}_j = \tilde{z}_{j-1} + \beta_{j-1} \tilde{p}_{j-1}. \qquad (12)$$

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

$$\alpha_j = \frac{\rho_{j-1}}{(\tilde{p}_j^T, q_j)}. \qquad (13)$$

Finally, one advances the source function according to

$$x_j = x_{j-1} + \alpha_j p_j, \qquad (14)$$

and the two residuals:

$$r_j = r_{j-1} - \alpha_j q_j, \qquad (15)$$

and,

$$\tilde{r}_j = \tilde{r}_{j-1} - \alpha_j \tilde{q}_j. \qquad (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

$$\frac{\|r\|_2}{\|b\|_2} < 10^{-2}. \qquad (17)$$

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

## 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_{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_{\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 $\epsilon$ values ranging from $10^{-4}$ to $10^{-12}$. Both the $\sqrt{\epsilon}$ surface law and thermalization lengths $\approx 1/\sqrt{\epsilon}$ are perfectly recovered, even for the numerically difficult cases of large albedos ($a \lesssim 1$).

The “true” error, that is,

$$T_e = \max \left( \frac{\left| S^{(j)} - S_{Edd} \right|}{S_{Edd}} \right), \qquad (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 $\epsilon = 10^{-12}$ in that case.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Solutions for the monochromatic scattering in an isothermal atmosphere, obtained for different values of $\epsilon$, ranging from $10^{-4}$ to $10^{-12}$. Both the expected $\sqrt{\epsilon}$ surface law and thermalization lengths $\approx 1/\sqrt{\epsilon}$ are perfectly recovered.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 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.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** History of the true error vs. iteration number for $\epsilon = 10^{-12}$, respectively, with Gauss-Seidel, SOR (with $\omega = 1.5$), and BiCG-P with the diagonal of A as a preconditioner.</figcaption></figure>

<table><caption>**Table 1.** Timing properties of the BiCG-P iterative scheme.</caption><tr><td>N<sub>τ</sub></td><td>5</td><td>8</td><td>11</td><td>14</td><td>17</td><td>20</td></tr><tr><td>t/t<sub>ALI</sub></td><td>0.85</td><td>1.2</td><td>1.8</td><td>2.1</td><td>2.3</td><td>2.7</td></tr></table>

## 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_τ$, for a $τ_{max} = 2 \times 10^8$, and $\epsilon = 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.

## 3.3. Sensitivity to the spatial refinement

We adopted the same slab model as the one reported in the previous section, with $\epsilon = 10^{-8}$, and we allowed the grid refinement to vary from 5 to 20 points per decade, with a step of 3. 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). 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^5$, sampled with a 9-point per decade grid, and $\epsilon = 10^{-6}$. The test consists in an 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
<figure><img src="image_6.png" /><figcaption>**Fig. 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.</figcaption></figure>

successful. In such a case, the timing properties of the *precondi-
tioned* 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 us-
ing the $A^T$ operator. A detailed trade-off analysis between over-
computing time and an improved convergence rate, with respect

to 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

Auer, L. H., & Paletou, F. 1994, A&A, 285, 675

Auer, L. H., Fabiani Bendicho, P., & Trujillo Bueno, J. 1994, A&A, 292, 599

Cannon, C. J. 1973, ApJ, 185, 621

Chevallier, L., Paletou, F., & Rutily, B. 2003, A&A, 411, 221

Cuny, Y. 1967, Ann. Ap., 30, 143

Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. H. 1997, A&A, 324, 161

Faurobert, M., Frisch, H., & Nagendra, K. N. 1997, A&A, 322, 896

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

Mihalas, D. 1978, Stellar Atmospheres (San Francisco: Freeman)

Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, JQSRT, 35, 431

Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171

Scharmer, G. B. 1981, ApJ, 249, 720

Saad, Y. 2008, Iterative methods for sparse linear systems (Philadelphia: SIAM)

Trujillo Bueno, J., & Fabiani Bendicho, P. 1995, ApJ, 455, 646

Trujillo Bueno, J., & Manso Sainz, R. 1999, ApJ, 516, 436

van Noort, M., Hubeny, I., & Lanz, T. 2002, ApJ, 568, 1066