A conjugate gradient method for the solution of the non-LTE line radiation transfer problem
aa r X i v : . [ a s t r o - ph . I M ] S e p Astronomy&Astrophysicsmanuscript no. arxivcg c (cid:13)
ESO 2018October 31, 2018
A conjugate gradient method for the solution of the non–LTE lineradiation transfer problem(Research Note)
F. Paletou and E. Anterrieu
Laboratoire d’Astrophysique de Toulouse–Tarbes, Universit´e de Toulouse, CNRS, 14 av. E. Belin, 31400 Toulouse, Francee-mail: [email protected]
Received May 14, 2009; accepted August 26, 2009
ABSTRACT
This study concerns the fast and accurate solution of the line radiation transfer problem, under non-LTE conditions. We proposeand evaluate an alternative iterative scheme to the classical ALI-Jacobi method, and to the more recently proposed Gauss-Seidel andSuccessive Over-Relaxation (GS / SOR) schemes. Our study is indeed based on the application of a preconditioned bi-conjugate gradi-ent method (BiCG-P). Standard tests, in 1D plane parallel geometry and in the frame of the two-level atom model, with monochromaticscattering, are discussed. Rates of convergence between the previously mentioned iterative schemes are compared, as well as theirrespective 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–LTEconditions, is a classical problem in astrophysics. Many numeri-cal methods have been used since the beginning of the computerera, from di ff erence 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-calledALI or Accelerated Λ -Iteration scheme is nowadays one ofthe most popular method for solving complex radiation transferproblems. Using as approximate operator the diagonal of the Λ operator (see e.g., Mihalas 1978), it is a Jacobi iterative scheme.It has been generalized for multilevel atom problems (Rybicki &Hummer 1991), multi-dimensional geometries (Auer & Paletou1994, Auer et al. 1994, van Noort et al. 2002), and polarizedradiation 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 non–LTE radiation transfer. They adaptedGauss-Seidel and Successive Over-Relaxation (GS / SOR) itera-tive schemes, well-known in applied mathematics as being su-perior , in terms of convergence rate, to the ALI–Jacobi iterativescheme. It is interesting to note that, besides their applicationto radiation transfer, both Jacobi and GS / SOR iterative schemeswere proposed during the xix th century.Conjugate gradient-like, hereafter CG-like, iterative methodswere proposed more recently by Hestenes & Stiefel (1952). Theyare of a distinct nature than ALI-Jacobi and GS / SOR methods.Unlike the later, CG-like methods are so-called non-stationary iterative methods (see e.g., Saad 2008). Very few articles dis-cussed 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 astrophysicalpurposes. Such an evaluation is the scope of the present research note.This alternative approach is relevant to the quest for ever fasterand accurate numerical methods for the solution of the non–LTEline radiation transfer problem.
2. The iterative scheme
In the two-level atom case , and with complete frequency redis-tribution, the non-LTE line source function is usually written as S ( τ ) = (1 − ε ) ¯ J ( τ ) + ε B ( τ ) , (1)where τ is the optical depth, ε is the collisional destruction prob-ability , B is the Planck function and ¯ J is the usual mean inten-sity:¯ J = I d Ω π Z + ∞−∞ φ ν I ν Ω d ν , (2)where the optical depth dependence of the specific intensity I ν Ω (and, possibly, of the line absorption profile φ ν ) is omitted forsimplicity.The mean intensity is usually written as the formal solutionof the radiative transfer equation¯ J = Λ [ S ] , (3)or, in other terms, the solution of the radiation transfer equationfor a known source function (see e.g., Mihalas 1978).Let us define A ≡ Id − (1 − ε ) Λ , the unknown x ≡ S and theright-hand side b ≡ ε B . Then, the solution of the radiative trans-fer equation is equivalent to solving the system of equations: Ax = b . (4) The albedo a = (1 − ε ) is more commonly used in general studiesof the radiation transfer equation. F. Paletou and E. Anterrieu: A conjugate gradient method for the solution of the non–LTE line radiation transfer problem (RN) Fig. 1.
Typical structure of the A operator, computed for the caseof a self-emitting 1D, plane parallel finite slab. Such an operatoris, more generally, non-symmetric positive definite. For the sakeof getting a better contrast for the o ff -diagonal parts, we magni-fied the negative values of the original operator by a factor of 5,and we display in fact the opposite of the modified operator.Among the various forms of CG-based methods, the bi-conjugate gradient (hereafter BiCG) method is more appropriatefor the case of operators A which are not symmetric positive def-inite. It is indeed the case for our radiative transfer problem, asillustrated in Fig. 1 where we display the structure of the oper-ator A for the case of a symmetric 1D grid of maximum opticaldepth τ max = × , discretized with 8 spatial points per decadeand using ε = − .A general description of the BiCG gradient method can befound in several classical textbooks of applied mathematics (seee.g., Saad 2008; it is important to note that this method requiresthe use of the transpose of the A operator, A T ). However, theBiCG method remains e ffi cient until the A operator becomesvery badly conditioned i.e., for very high albedo cases, typi-cally when ε < − . For the later cases, it is therefore crucial toswitch to a more e ffi cient scheme, using preconditioning of thesystem of equations.Hereafter, we derive the algorithm for the BiCG method withpreconditioning (hereafter BiCG-P). In such a case, one seeksinstead for a solution of: M − Ax = M − b , (5)where the matrix M should be easy to invert. A “natural” choicefor the non-LTE radiation transfer problem, is to preconditionthe system of equations with the diagonal of the full operator A .From an initial guess for the source function, x = S (0) , com-pute a residual r such as: r = b − Ax . (6)Set the second residual ˜ r such that ( r , ˜ r ) ,
0, where ( a , b )means the inner product between vectors a and b . For all thecases considered by us, setting r = ˜ r proved to be adequate.Now, the BiCG-P algorithm consists in running until con-vergence the following iterative scheme, where j is the iterativescheme index. The first steps consist in solving: Mz j − = r j − , (7) and M T ˜ z j − = ˜ r j − . (8)These steps are both simple and fast to compute, since M is adiagonal operator. Then one computes: ρ j − = ( z Tj − , ˜ r j − ) . (9)If ρ j − = j = p j = z j − and ˜ p j = ˜ z j − . For j >
1, compute: β j − = ρ j − /ρ j − , (10) p j = z j − + β j − p j − , (11)and:˜ p j = ˜ z j − + β j − ˜ p j − . (12)Then make q j = Ap j and ˜ q j = A T ˜ p j , as well as: α j = ρ j − ( ˜ p Tj , q j ) . (13)Finally, one advances the source function according to: x j = x j − + α j p j , (14)and the two residuals: r j = r j − − α j q j , (15)and,˜ r j = ˜ r j − − α j ˜ q j . (16)This process, from Eq. (7) to Eq. (16), is then repeated toconvergence. We have found that a convenient stopping criterionis to interrupt the iterative process when: k r k k b k < − . (17)Indeed, this stopping criterion guarantees that the floor value ofthe true error, defined in §
3. Results
The solution of the non–LTE radiative transfer problem for aconstant property, plane-parallel 1D slab is a standard test forany new numerical method. Indeed, the so-called monochro-matic scattering case allows a direct comparison of the numer-ical solutions to the analytical “solution”, S Edd , which can bederived in Eddington’s approximation (see e.g., the discussionin §
5. of Chevallier et al. 2003). . Paletou and E. Anterrieu: A conjugate gradient method for the solution of the non–LTE line radiation transfer problem (RN) Fig. 2.
Solutions for the monochromatic scattering in an isother-mal atmosphere, obtained for di ff erent values of ε , ranging from10 − to 10 − . Both the expected √ ε surface law and thermali-sation lengths ≈ / √ ε are perfectly recovered. In Fig. 2, we display the run of the source function with opticaldepth, for a self-emitting slab of total optical thickness τ max = × , using a 9-point per decade spatial grid, a one-point angularquadrature with µ = ± / √
3, and ε values ranging from 10 − to10 − .Both the √ ε surface law and thermalisation lengths ≈ / √ ε are perfectly recovered, even for the numerically di ffi cult casesof large ( a /
1) albedos.The “true” error, that is: T e = max | S ( j ) − S Edd | S Edd ! , (18)is displayed in Fig. 3, respectively for the GS, SOR and BiCG-Piterative schemes. SOR was used with a ω = . ε = − in that case. The typical timing properties of the BiCG-P scheme, for thecases considered here, are summarized in Table 1. While eachGS / SOR iteration is about 20% longer than one ALI-iteration,the increase of time per iteration, with the number of depthpoints, for BiCG-P is significant. However, the balance betweencomputing time and reduction of the number of iterations alwaysremains in favour of the BiCG-P scheme vs. GS / SOR and, a for-tiori, ALI.
We adopted the same slab model as the one reported in the pre-vious section, with ε = − , and we allowed the grid refinementto vary from 5 to 20 points per decade, with a step of 3.Figure 4 shows the sensitivity of the method with grid re-finement. Increasing grid refinement corresponds to a decreasingvalue of the floor value of the true error, which demonstrates the Fig. 3.
History of the true error vs. iteration number for ε = − , respectively with Gauss-Seidel, SOR (with ω = .
5) andBiCG-P with the diagonal of A as a preconditioner.
Table 1.
Timing properties of the BiCG-P iterative schemes, de-pending on the optical depth grid refinement expressed as thenumber of point per decade, N τ , for a τ max = × , and ε = − self-emitting slab. We use a time normalized to thetime for a ALI-iteration. Note that GS / SOR is, in comparison,about 20% slower than ALI. N τ / t ALI need of a su ffi cient spatial resolution for the sake of accuracy onthe numerical solutions. The number of iterative steps necessaryto fulfill the stopping criterion defined in Eq. (17) is practicallylinear in the number of depth points. A strong argument in favour of the Gauss-Seidel iterativescheme resides in its smoothing capability. This point is partic-ularly important in the frame of multi-grid methods for the so-lution of non–LTE radiative transfer problems (see e.g., FabianiBendicho et al. 1997).We have reproduced the test initially proposed by TrujilloBueno & Fabiani Bendicho (1995, see their Fig. 9) for the GSmethod. Here we consider a semi-infinite slab of maximum op-tical thickness 10 , sampled with a 9-point per decade grid, and ε = − . The test consists in injecting ab initio a high-frequencycomponent into the error, S (0) − S ( ∞ ) , on the initial guess for thesource function. In Fig. 5, the initial error we used is easily rec-ognized as the hyperbolic tangent-like curve, to which a high-frequency component was added.After the very first iterate of the BiCG-P iterative process,the error has an amplitude of ≈ .
75 and, more importantly, itis already completely smoothed. Other error curves displayed inFig. 5, of continuously decreasing amplitudes, are the successiveerrors at iterates 2 to 5.This simple numerical experiment demonstrates that BiCG-Pcan also be used as an e ffi cient smoother for multi-grid methods. F. Paletou and E. Anterrieu: A conjugate gradient method for the solution of the non–LTE line radiation transfer problem (RN)
Fig. 4.
Evolution of the true error for varying spatial grid refine-ments, from 5 to 20 points per decade, with a step of 3. Thevalue of the true error plateau decreases with an increasing gridrefinement.
Fig. 5.
Initialization and evolution of the error, S ( j ) − S ( ∞ ) , duringthe 5 first iterations of the BiCG-P scheme. The high-frequencycomponent introduced ab initio in the source function is re-moved after the very first iteration of BiCG-P. Other curves, ofcontinuously decreasing amplitudes, correspond to iterations 2to 5.
4. Conclusion
We propose an alternative method for the solution of the non–LTE radiation transfer problem. Preliminary but standards testsfor the two-level atom case in 1D plane parallel geometry aresuccessful. In such a case, the timing properties of the precondi-tioned
BiCG method are comparable to the ones of ALI, GS andSOR stationnary iterative methods. However the convergencerate of BiCG-P is signifiantly better.The main potential drawback of the BiCG-P method is inthe usage of the A T operator. A detailed trade-o ff analysis be- tween such over-computing time and an improved convergencerate, with respect to the ones of GS / SOR iterative schemes, iscurrently being conducted.Multi-level and multi-dimensional cases will be consideredfurther. However, the use of BiCG-P does not require the cum-bersome modifications of multi-dimensional formal solvers re-quired by GS / SOR methods (see e.g., L´eger et al. 2007).
Acknowledgements.
We are grateful to Drs. Bernard Rutily, and Lo¨ıc Chevallierfor numerous discussions about the solution of the radiative transfer equation.