Issue 
A&A
Volume 498, Number 3, May II 2009



Page(s)  981  985  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/200911661  
Published online  12 March 2009 
A 3D radiative transfer framework
IV. Spherical and cylindrical coordinate systems
P. H. Hauschildt^{1}  E. Baron^{1,2,3}
1 
Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg, Germany
2 
Dept. of Physics and Astronomy, University of
Oklahoma, 440 W. Brooks, Rm 100, Norman, OK 73019, USA
3 
Computational Research Division, Lawrence Berkeley National Laboratory,
MS
50F1650, 1 Cyclotron Rd, Berkeley, CA 947208139, USA
Received 14 January 2009 / Accepted 8 March 2009
Abstract
Aims. We extend our framework for 3D radiative transfer calculations with a nonlocal operator splitting methods along (full) characteristics to spherical and cylindrical coordinate systems. These coordinate systems are better suited to a number of physical problems than Cartesian coordinates.
Methods. The scattering problem for line transfer is solved via means of an operator splitting (OS) technique. The formal solution is based on a full characteristics method. The approximate
operator is constructed considering nearest neighbors exactly. The code is parallelized over both wavelength and solid angle using the MPI library.
Results. We present the results of several test cases with different values of the thermalization parameter for the different coordinate systems. The results are directly compared to 1D plane parallel tests. The 3D results agree very well with the welltested 1D calculations.
Conclusions. Advances in modern computers will make realistic 3D radiative transfer calculations possible in the near future.
Key words: radiative transfer  stars: atmospheres
1 Introduction
In Hauschildt & Baron (2006, hereafter: Paper I), Baron & Hauschildt (2007, hereafter: Paper II) and (Hauschildt & Baron 2008, hereafter: Paper III) we described a framework for the solution of the radiative transfer equation for scattering continua and lines in 3D (when we say 3D we mean three spatial dimensions, plus three momentum dimensions) for the time independent, static case using full characteristics based on voxels in Cartesian coordinates. Cartesian coordinates are by no means an optimal choice for many 3D physical systems, for example supernova atmospheres are better described in 3D spherical and accretion disks are better described in 3D cylindrical coordinates. These coordinate systems make it a little harder and computationally more expensive to track the characteristics through the curvilinear voxel grid. However, this is more than overcompensated by the higher efficiency of voxel usage, for example consider putting a roughly spherical structure where density and composition vary strongly with radius into a Cartesian coordinate system cube.
We describe our method, its rate of convergence, and present comparisons to our welltested 1D calculations.
2 Method
In the following discussion we use notation of Papers IIII. The basic framework and the methods used for the formal solution and the solution of the scattering problem via operator splitting are discussed in detail in Papers I and II and will not be repeated here.
We use a ``full characteristics'' approach that completely tracks a set of characteristics of the radiative transfer equation from the inner or outer boundary through the computational domain to their exit voxels and takes care that each voxel is hit by at least one characteristic per solid angle. One characteristic is started at each boundary voxel and then tracked through the voxel grid until it leaves the other boundary. The direction of a bundle of characteristics is determined by a set of solid angles which correspond to a normalized momentum space vector (p_{x},p_{y},p_{z}).
In contrast to the case of Cartesian coordinates the volumes of the voxels can be vastly different in spherical and Cartesian coordinate systems. Therefore, tracking a characteristic across the voxel grid is done using an adaptive stepsize control along the characteristic so that a characteristic steps from one voxel to its neighbor in the direction of the characteristic. In addition, we make sure that each voxel is hit by at least one characteristic for each solid angle, which is important for the inner (small radii) voxels.
The code is parallelized as described in Papers I and II.
3 Spherical coordinates
3.1 Testing environment
We use the framework discussed in Papers I and II as the baseline for the line transfer problems discussed in this paper. Our basic setup is similar to that discussed in Paper II. We use a sphere with a grey continuum opacity parameterized by a power law in the continuum optical depth . The basic model parameters are 1.
 Inner radius cm, outer radius cm (large test case) or cm (small test case).
 2.
 Minimum optical depth in the continuum and maximum optical depth in the continuum .
 3.
 Constant temperature with K, therefore, the Planck function B is also constant.
 4.
 Outer boundary condition and diffusion inner boundary condition for all wavelengths.
 5.
 Continuum extinction , with the constant Cfixed by the radius and optical depth grids, for the continuum tests, for the line transfer tests
 6.
 Parameterized coherent & isotropic continuum scattering by
defining
with . and are the continuum absorption and scattering coefficients.
The line of the simple 2level model atom is parameterized by the ratio of the profile averaged line opacity to the continuum opacity and the line thermalization parameter . For the test cases presented below, we have used , , and a grey temperature structure. We set to simulate a strong line, with varying (see below). With this setup, the optical depths as seen in the line range from 10^{2} to 10^{6}. We use 32 wavelength points to model the full line profile, including wavelengths outside the line for the continuum. We did not require the line to thermalize at the center of the test configurations, this is a typical situation for a full 3D configuration as the location (or even existence) of the thermalization depths becomes more ambiguous than in the 1D case.
In the following we discuss the results of various tests. Unless otherwise stated, the tests were run on parallel computers using 4 to 512 sCPUs (depending on the test). In spherical coordinates, we use a voxel grid with voxels; this corresponds to an angular resolution of about 5.5 degrees in and . Projected on the surface of the Earth, this means that Europe is covered by about 43 voxels. The spherical test cases use incoming intensities of zero at the outside ( and I=B at the inside ( ).
Figure 1: The ratio of the mean intensity J and the Planck function B as function of the radial optical depth for the ``large'' test case in spherical coordinates and different values . The optical depth is measured along the radius r. 

Open with DEXTER 
3.2 Results
3.2.1 Continuum tests
In Fig. 1 we show the results for 3D continuum transfer with 3D spherical coordinates for various values of the thermalization parameter . The results are compared to 1D results (full line). In all cases, the agreement is excellent, even for the case with . Compared to Paper I the agreement is better, this is caused by the better discretization of a spherical configuration in spherical coordinates, rather than fitting a sphere into a Cartesian grid. It illustrates the importance of adapted coordinate systems to obtain the optimal mix between computational efficiency and numerical accuracy.
Good coverage of the solid angle is important for 3D transfer. The 1D code automatically adjusts the angle coverage related to the geometry of the configuration. However, in the 3D code the solid angles have to be prescribed. In Fig. 2 we demonstrate the importance of solid angle sampling for a ``small'' test case with . The plots show the ratio J/B as a function of radial optical depth for various solid angle samplings and for all compared to the 1D results. It is immediately clear that for the ``small'' test case a solid angle sampling of 64^{2} points or more is reasonable, but that 16^{2} or less points result in a considerable scatter of the results. This result does not depend strongly on the radial extension of the configuration, as shown in Fig. 3.
The convergence rate of the nearest neighbor (3^{3} neighbor voxels are considered in the , see Paper I for details) with Ng acceleration is virtually identical to that of the 1D tridiagonal case, cf. Fig. 4. The importance of a nonlocal for good convergence in the 3D case can be seen by comparing the results obtained with a local (diagonal) . Its convergence with Ng acceleration is actually worse than that of the nonlocal without Ng acceleration. The excellent convergence of the nonlocal can be used to reduce the computing time by using a less stringent convergence criterion that still gives results accurate enough for a particular application. One has to be careful to be not too aggressive in relaxing the convergence criterion. This is highlighted in Fig. 5 where the results for the large test case are shown for different combinations of with Ng acceleration and difference convergence limits. It is obvious that the combination of Ng acceleration and a relatively large convergence limit with a local (diagonal) is not acceptable, whereas the same convergence limit (10^{2}maximum relative change of J) with a nonlocal (nearest neighbor) gives usable results even without Ng acceleration. Setting a convergence limit of 10^{4} gives reasonable results in all test cases and can reduce computing time by about 50%.
Figure 2: Effect of solid angle resolution on the results of the 3D calculations. The plots show the results for the ``small'' continuum test case in spherical coordinates with and for various solid angle resolutions with from (256,256) to (8,8). Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER 
Figure 3: Effect of solid angle resolution on the results of the 3D calculations in spherical coordinates. The plots show the results for the ``large'' continuum test case with and for various solid angle resolutions with from (256,256) to (8,8). Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER 
Figure 4: Convergence behavior for the ``large'' continuum test case in spherical coordinates with for various combinations of local/nonlocal and Ng acceleration. ``NN'' indicates the nonlocal nearest neighbor operator, ``diagonal'' the local . For comparison the results for the 1D spherical code and the 3D lambda iterations are also plotted. 

Open with DEXTER 
Figure 5: Effects of relaxing the convergence criterions on the results of the 3D calculations in spherical coordinates. The plots show the results for the ``large'' continuum test case with and for various convergence limits for nearest neighbor (NN) and diagonal (Diag) operators with and without Ng acceleration. The numbers are the maximum allowed relative change over all voxels used to halt the iterations. Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER 
Figure 6: The ratio of the line mean intensity and the Planck function B as a function of the radial optical depth for the ``large'' test case in spherical coordinates for different values of . The optical depth is measured along the radius r. 

Open with DEXTER 
Figure 7: Flux vectors in the outermost voxels for the line transfer test in spherical coordinates with and 64^{2} solid angle points. The top panel shows the radial component F_{r} of the voxel flux vector compared to the results for the 1D comparison calculation. The middle and bottom panels show the ratio of the polar and azimuthal flux components to F_{r}, respectively. 

Open with DEXTER 
Figure 8: Flux vectors in the outermost voxels for the line transfer test in spherical coordinates with and 512^{2} solid angle points. The top panel shows the radial component F_{r} of the voxel flux vector compared to the results for the 1D comparison calculation. The middle and bottom panels show the ratio of the polar and azimuthal flux components to F_{r}, respectively. 

Open with DEXTER 
With 128 MPI processes on the SGI ICE system of the Höchstleistungs Rechenzentrum Nord (HLRN) the continuum test with and and spatial points and 256^{2} solid angle points requires 16 iterations and 1543 s wallclock time, where a formal solution takes about 80 s, the construction of the nonlocal about 140 s and about 1.6 s are used for the solution of the operator splitting step. For the same case with 64^{2} solid angles the total wallclock time is 230 s. The memory requirements are about 45 MB/process.
3.2.2 Line tests
In Fig. 6 we show the results of a comparison of the ratio of the line mean intensities and the Planck function Bfor the ``large'' test case and various values of . The 1D comparison models assume a diffusive inner boundary condition, whereas the 3D models use I=B as the boundary condition for the inner boundary. Therefore, the models will differ close to the inner boundary. Away from the inner boundary the results of the 1D and 3D calculations agree very well with each other. The number of iterations required to reach a prescribed accuracy is similar to the continuum tests described above. This shows that the quality of the 3D calculations match the 1D results, which is very important for future applications.
In Figs. 7 and 8 we show the flux vectors of the outermost voxels (outermost radius, all points) for the line transfer case with and for different numbers of solid angle points. In the test case, the radial component of the flux vector F_{r} should be the only nonzero component as in spherical symmetry both and are zero. The figures show that this is indeed the case if a large enough number of solid angle points is used in the calculation. If the number is too small, and can have errors of about 10% and F_{r} scatters considerably around the 1D comparison. For 512^{2} solid angle points, the scatter of F_{r} is too small to be seen in the plot and and are smaller than 1% of F_{r} over the whole surface of the spherical grid. This highlights the necessity of having enough solid angle points to obtain good accuracy. Figure 9 shows F_{r} compared to the 1D results for a number of values of . In all cases the results agree very well with each other, showing the good numerical accuracy of the 3D code. The results of these tests are valid only for the test cases. Therefore, a similar test should be performed for each ``real'' 3D configuration individually in order to make sure that enough solid angles are used for the accuracy requirements of each individual problem.
Figure 9: Radial components F_{r} of the flux vectors in the outermost voxels for the line transfer tests in spherical coordinates with and 256^{2} solid angle points. The results of the 1D calculations are shown for comparison. 

Open with DEXTER 
With 512 MPI processes on the SGI ICE system of the HLRN the line test with and and spatial points, 32 wavelength points and 256^{2} solid angle points requires 30 iterations and 20 287 s wallclock time and for the same case with 64^{2}solid angles the total wallclock time is 1400 s. In both cases we used 32 ``clusters'' with 16 MPI processes each so that for each wavelength 16 MPI processes work on a formal solution (see Paper II for details).
4 Cylindrical coordinates
The setup for the test of the code in cylindrical coordinates is essentially identical to the tests of the spherical coordinate system version. However, describing a spherically symmetric objects in cylindrical coordinates is nonoptimal and we expect results comparable to those presented in Papers I and II. The cylindrical coordinate system version of the code allows for different maps in and zfor flexibility describing objects with nearly cylindrical symmetry.
The differences in the basic setup compared to the spherical coordinates tests are
 1.
 inner radius cm, outer radius cm;
 2.
 linear radial grid;
 3.
 Continuum extinction , with the constant Cfixed by the radius and optical depth grids, for the continuum tests, for the line transfer tests.
4.1 Results
For the cylindrical coordinate tests we show just the results of a standard test discussed in Paper I and above for the spherical case with . The results for the continuum test case are presented in Fig. 10. The results agree reasonably well with the 1D test case, comparable to the results described in Paper I. The reason for this is that describing a sphere in cylindrical coordinates is not an ideal use of this coordinate system, therefore much higher grid resolution is required to reach reasonable accuracy (the results for grids with are substantially worse). As for the test cases in spherical coordinates, the resolution in solid angles is very important for accurate results. The convergence behavior is the same as in the case of spherical coordinates.
5 Conclusions
We have discussed the extension of our 3D framework to 3D spherical and cylindrical coordinate systems. These coordinates are better suited than Cartesian coordinates for problems that are approximately spherical (e.g., supernova envelopes) or cylindrical (e.g., protoplanetary disks). The convergence properties and the accuracy of the solutions are comparable to the 1D solution for the test cases discussed here. These additional coordinate systems allow for more general use with higher internal accuracy while using less memory than the purely Cartesian coordinate system.
Figure 10: Results for an example of a grey radiative transfer calculation for a sphere in cylindrical 3D coordinates for with and 64^{2} solid angle points. The results of the 1D calculations are shown for comparison with * symbols. 

Open with DEXTER 
The next step in the development of the 3D framework is to add the treatment of velocity fields in the observer's frame and the comoving frame. The latter will be important for applications with relativistic velocity fields, e.g., supernovae, while the former is interesting for modeling the spectra of convective flow models.
Acknowledgements
This work was supported in part by by NASA grants NAG53505 and NAG512127, NSF grants AST0307323, and AST0707704, and US DOE Grant DEFG0207ER41517, as well as DFG GrK 1351 and SFB 676. Some of the calculations presented here were performed at the Höchstleistungs Rechenzentrum Nord (HLRN); at the NASA's Advanced Supercomputing Division's Project Columbia, at the Hamburger Sternwarte Apple G5 and Delta Opteron clusters financially supported by the DFG and the State of Hamburg; and at the National Energy Research Supercomputer Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC0376SF00098. We thank all these institutions for a generous allocation of computer time.
References
 Baron, E., & Hauschildt, P. H. 2007, A&A, 468, 255 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Hauschildt, P. H. 1992, JQSRT, 47, 433 [NASA ADS] (In the text)
 Hauschildt, P. H., & Baron, E. 2006, A&A, 451, 273 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
 Hauschildt, P. H., & Baron, E. 2008, A&A, 490, 873 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
All Figures
Figure 1: The ratio of the mean intensity J and the Planck function B as function of the radial optical depth for the ``large'' test case in spherical coordinates and different values . The optical depth is measured along the radius r. 

Open with DEXTER  
In the text 
Figure 2: Effect of solid angle resolution on the results of the 3D calculations. The plots show the results for the ``small'' continuum test case in spherical coordinates with and for various solid angle resolutions with from (256,256) to (8,8). Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER  
In the text 
Figure 3: Effect of solid angle resolution on the results of the 3D calculations in spherical coordinates. The plots show the results for the ``large'' continuum test case with and for various solid angle resolutions with from (256,256) to (8,8). Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER  
In the text 
Figure 4: Convergence behavior for the ``large'' continuum test case in spherical coordinates with for various combinations of local/nonlocal and Ng acceleration. ``NN'' indicates the nonlocal nearest neighbor operator, ``diagonal'' the local . For comparison the results for the 1D spherical code and the 3D lambda iterations are also plotted. 

Open with DEXTER  
In the text 
Figure 5: Effects of relaxing the convergence criterions on the results of the 3D calculations in spherical coordinates. The plots show the results for the ``large'' continuum test case with and for various convergence limits for nearest neighbor (NN) and diagonal (Diag) operators with and without Ng acceleration. The numbers are the maximum allowed relative change over all voxels used to halt the iterations. Full line: results of the 1D calculation, dots: results from the 3D calculations. 

Open with DEXTER  
In the text 
Figure 6: The ratio of the line mean intensity and the Planck function B as a function of the radial optical depth for the ``large'' test case in spherical coordinates for different values of . The optical depth is measured along the radius r. 

Open with DEXTER  
In the text 
Figure 7: Flux vectors in the outermost voxels for the line transfer test in spherical coordinates with and 64^{2} solid angle points. The top panel shows the radial component F_{r} of the voxel flux vector compared to the results for the 1D comparison calculation. The middle and bottom panels show the ratio of the polar and azimuthal flux components to F_{r}, respectively. 

Open with DEXTER  
In the text 
Figure 8: Flux vectors in the outermost voxels for the line transfer test in spherical coordinates with and 512^{2} solid angle points. The top panel shows the radial component F_{r} of the voxel flux vector compared to the results for the 1D comparison calculation. The middle and bottom panels show the ratio of the polar and azimuthal flux components to F_{r}, respectively. 

Open with DEXTER  
In the text 
Figure 9: Radial components F_{r} of the flux vectors in the outermost voxels for the line transfer tests in spherical coordinates with and 256^{2} solid angle points. The results of the 1D calculations are shown for comparison. 

Open with DEXTER  
In the text 
Figure 10: Results for an example of a grey radiative transfer calculation for a sphere in cylindrical 3D coordinates for with and 64^{2} solid angle points. The results of the 1D calculations are shown for comparison with * symbols. 

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.