$\epsilon =10^{-4}$

For this test, we use the same basic structure setup as for the LTE test, however we now use $\epsilon =10^{-4}$ for a scattering dominated atmosphere. The comparison with the results obtained with the 1D RT code (with 100 radial points) are compared to the 3D RT code results in Fig. 3. Note the factor of about 1000 difference between the solution for $\epsilon =10^{-4}$ and the results of the formal solution with $S=B$ shown in Fig. 2 at the largest distances (the iterations were started with $S=B$). Figure 3 shows the results for the `large' and `medium' spatial grids and for 2 different solid angle discretizations ($64^2$ and $16^2$ solid angle points). For both the 1D and the 3D code the mean intensities were iterated to a relative accuracy of $10^{-8}$ (see below for details on convergence rates). The graph highlights the need for a rather fine solid angle grid for the test problem. With a small number of angular points (bottom panel in Fig. 3) the numerically induced spread of the mean intensities is significantly larger than with a $4^2$ finer angular grid (shown in the middle panel). A change in spatial grid resolution has a much smaller effect on the results than the number of solid angle points as shown in the top two panels of Fig. 3. The spatial resolution of the `small' grid is, however, too coarse to represent the test problem, in particular in the inner parts of the structure.

The importance of the angular resolution is also demonstrated in Fig. 4 which shows the contours of the mean intensity on the six `faces' of the voxel cube for $129^3$ spatial points and $16^2$ (left panel) and $64^2$ (right panel) solid angles. The calculation with $64^2$ solid angle points shows symmetric contours whereas the smaller $16^2$ angle points model shows asymmetries and `banding' like structures on all faces (in particular on the $x-y$ faces). This clearly indicates that for the test conditions used the angular resolution has to be larger than $16^2$ for accurate solutions. This can also be seen in surface plots of the mean intensities at the $-n_{\rm x}$ faces of the $z-y$ planes for both calculations, cf. Fig. 5. The surface calculated with $64^2$ angles is much smoother and shows no or little banding compared to the surface for $16^2$ angular points. The effects of higher spatial resolution on the $J$ surface at the $z-y$ face is shown in Fig. 6. Clearly, $64^2$ angles produce a smooth surface without significant artifacts for the test case.

The convergence properties of the various methods for this test case are shown in Fig. 7. As expected, the $\Lambda$ iteration converges very slowly and requires more than 300 iterations to reach a relative error of less than $10^{-8}$. The diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operator converges significantly better than the $\Lambda$ iteration but is still rather slowly convergent. The nearest neighbor $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ converges substantially faster (by more than a factor of 3 than the diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$). The diagonal and nearest neighbor iterations can be accelerated by using Ng's method (Ng, 1974), as shown in Fig. 7. In both cases, a 4th order Ng acceleration was used after 20 initial iterations, starting Ng acceleration too early can lead to convergence failures since Ng acceleration is based on extrapolation. An attempt to apply the Ng acceleration to the $\Lambda$ iteration was not successful, similar to the 1D case, as the convergence rate of the $\Lambda$ iteration appears to be too small to be useful for the Ng method. For comparison, we show in Fig. 7 the convergence properties of the 1D RT solver with a tri-diagonal $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ and Ng acceleration (here started much earlier). The convergence rates of the 1D tri-diagonal and 3D nearest neighbor methods are very comparable. The convergence properties are also relevant for the overall speed of the method: whereas the time for the formal solution (for a given number of voxels and processors) depends roughly linearly on the total number of solid angle points, the time for the solution of the operator splitting step does not depend on the number of angle points and actually decreases with iteration number if the Jordan, Gauss-Seidel or Rapido solvers are used. Therefore, the nearest neighbor operator will become more and more efficient as the number of angles and/or voxels increases.

The smaller initial corrections of the $\Lambda$ iteration are actually an indication that it just corrects too little, the operator splitting method gives initially much larger corrections. This is exactly like similar test cases in 1D that we have tried. The initial corrections are so large because the initial guesses are very wrong (intentionally) so that between initial guess and final result we have changes by close to 10 orders of magnitude. This means that the initial corrections of the operator splitting method have to be large.

The convergence rates will depend on the grid resolution as the optical depths between voxels will be smaller for larger resolutions and thus the coupling between voxels will be stronger. This effect is shown in Fig. 8 where we show the convergence rates for various grids. The convergence rates are independent of the number of solid angle points but depend weakly on the number of grid points, as expected. Thus, for voxel grids with higher resolution, a larger $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ (more neighbors) will be more useful than for coarser voxel grids.

Figure 9 shows the convergence rates for the $\epsilon =10^{-4}$ test case for different $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ bandwidths. The wider $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$'s lead to improved convergence rates; however, the net wallclock time increases for the wider $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$'s in the parallel code. This is caused by the substantially larger amount of information that needs to be passed between processes in order to build the wider $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$. In addition, the memory requirements for a given number of voxels scales like the cube of the number of neighbors considered in the $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$. As Fig. 9 demonstrates, the convergence is strongly improved by the use of Ng acceleration, however, for larger $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$'s the Ng method can be detrimental for wider $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ compared to narrower $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$.