$\epsilon =10^{-8}$

The final test we present in this paper considers a case with a much larger scattering contribution: $\epsilon =10^{-8}$. The results for a grid with $129^3$ voxels and $64^2$ angular points are shown in Fig. 10. Fig. 11 shows the results for slices along the coordinate axes (the two coordinates being centered, respectively). Even though the dynamic range of the mean intensities is huge, nearly 12 dex, the results are quite accurate. For this case the lack of spatial resolution in the inner parts of the voxel grid can be seen. Here $\vert\nabla J\vert$ is huge and cannot be fully resolved by the 3D RT code (the 1D code naturally has much higher resolution). However, only a few voxels away from the center the results agree very well.

The convergence plots in Fig 12 show the results for a very difficult test case with $\tau _{\rm std}^{\rm max} = 10^{8}$ for a fixed voxel and solid angle grid with $65^3$ voxels and $16^2$ angular points. For this test, the $\Lambda$ iteration fails completely. The diagonal operator provides significant speed-up, but still requires more than 1600 iterations to reach the required convergence limit. For this test, the Ng acceleration does not work with the diagonal operator, the iteration process failed immediately after it was started. It is likely that a better result could be obtained if Ng acceleration is started in the steep part of the diagonal operator's convergence, e.g., after about 500 regular iterations. The nearest neighbor operator leads to much faster convergence, even without Ng acceleration the solution converges in about 450 iterations. Here, Ng acceleration works very well with the nearest neighbor operator, convergence is reached after 177 iterations. This is still about a factor of 2 more than for the 1D code, but much better than with the diagonal operator. The variation of the convergence rate for the nearest neighbor operator and Ng acceleration with the size of the solid angle grid is shown in Fig. 13. The case with the smallest angular grid ($16^2$ points) actually converges more slowly than the $32^2$ and $64^2$ grids. The higher resolution grids convergence rate compares well with the 1D code. This highlights the importance of the non-local $\ifmmode{\Lambda^*}\else\hbox{$\Lambda^*$}\fi$ operator and a large enough solid angle grid for rapid convergence and accuracy.