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/0004-6361/200911661 | |
Published online | 12 March 2009 |
A 3D radiative transfer framework
IV. Spherical and cylindrical coordinate systems
P. H. Hauschildt1 - E. Baron1,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
50F-1650, 1 Cyclotron Rd, Berkeley, CA 94720-8139, USA
Received 14 January 2009 / Accepted 8 March 2009
Abstract
Aims. We extend our framework for 3D radiative transfer calculations with a non-local 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 well-tested 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 well-tested 1-D calculations.
2 Method
In the following discussion we use notation of Papers I-III. 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
(px,py,pz).
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
- 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 2-level 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 106. 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
|
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
642 points or more is reasonable, but that 162 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
(33 neighbor voxels
are considered in the
,
see Paper I for details) with Ng acceleration
is virtually identical to that of the 1D tri-diagonal
case, cf.
Fig. 4. The importance of a non-local
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 non-local
without Ng acceleration.
The excellent convergence of the non-local
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-2maximum relative change of J) with a non-local (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
|
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
|
Open with DEXTER |
![]() |
Figure 4:
Convergence behavior for the
``large'' continuum
test case in spherical coordinates with
|
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
|
Open with DEXTER |
![]() |
Figure 6:
The ratio of the line mean intensity |
Open with DEXTER |
![]() |
Figure 7:
Flux vectors in the outermost voxels
for the line transfer test in spherical coordinates with
|
Open with DEXTER |
![]() |
Figure 8:
Flux vectors in the outermost voxels
for the line transfer test in spherical coordinates with
|
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
2562 solid angle points requires 16 iterations and 1543 s wallclock time,
where a formal solution takes about 80 s, the construction of the non-local
about 140 s and about 1.6 s are used for the solution of the
operator splitting step. For the same case with 642 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 Fr should be the only non-zero 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 Fr scatters considerably around the 1D comparison.
For 5122 solid angle points, the scatter of Fr is too small to be seen in
the plot and
and
are smaller than 1% of Fr 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 Fr 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 Fr of the flux vectors in the outermost voxels
for the line transfer tests in spherical coordinates with
|
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 2562 solid angle points requires
30 iterations and 20 287 s wallclock time and for the same case with 642solid 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 non-optimal 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
|
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 co-moving 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 NAG5-3505 and NAG5-12127, NSF grants AST-0307323, and AST-0707704, and US DOE Grant DE-FG02-07ER41517, 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. DE-AC03-76SF00098. 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
|
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
|
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
|
Open with DEXTER | |
In the text |
![]() |
Figure 4:
Convergence behavior for the
``large'' continuum
test case in spherical coordinates with
|
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
|
Open with DEXTER | |
In the text |
![]() |
Figure 6:
The ratio of the line mean intensity |
Open with DEXTER | |
In the text |
![]() |
Figure 7:
Flux vectors in the outermost voxels
for the line transfer test in spherical coordinates with
|
Open with DEXTER | |
In the text |
![]() |
Figure 8:
Flux vectors in the outermost voxels
for the line transfer test in spherical coordinates with
|
Open with DEXTER | |
In the text |
![]() |
Figure 9:
Radial components Fr of the flux vectors in the outermost voxels
for the line transfer tests in spherical coordinates with
|
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
|
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.